**Flow–Sediment Turbulent Ejections: Interaction between Surface and Subsurface Flow in Gravel-Bed Contaminated by Fine Sediment**

**Bustamante-Penagos N. 1,\* and Niño Y. 1,2**


Received: 3 March 2020; Accepted: 23 April 2020; Published: 3 June 2020

**Abstract:** Several researchers have studied turbulent structures, such as ejections, sweeps, and outwards and inwards interactions in flumes, where the streamwise velocity dominates over vertical and transversal velocities. However, this research presents an experimental study in which there are ejections associated with the interchange between surface and subsurface water, where the vertical velocity dominates over the streamwise component. The experiment is related to a surface alluvial stream that is polluted with fine sediment, which is percolated into the bed. The subsurface flow is modified by a lower permeability associated with the fine sediment and emerges to the surface current. Quasi-steady ejections are produced that drag fine sediment into the surface flow. Particle image velocimetry (PIV) measured the velocity field before and after the ejection. The velocity data were analyzed by scatter plots, power spectra, and wavelet analysis of turbulent fluctuations, finding changes in the distribution of turbulence interactions with and without the presence of fine deposits. The flow sediment ejection changes the patterns of turbulent structures and the distribution of the turbulence interactions that have been reported in open channels without subsurface flows.

**Keywords:** ejections; turbulence interactions; gravel beds; sediment transport; surface and subsurface flows

#### **1. Introduction**

Landslides, volcanos, or anthropogenic changes may modify the availability of fine sediment in rivers, reservoirs, or lakes [1]. Fine sediment can cause pollution in gravel beds. This contamination has a high environmental impact because the porosity of the gravel beds is a reservoir for the deposition of fine sediment [2]. The intrusion of fine sediment into gravel bed streams generates changes in the hyporheic exchange, nutrients cycling, low oxygenation of fish eggs, etc. [3–5]. Additionally, the hyporreic zone also have an important coupling between the subsurface groundwater system and surface water, such as rivers or lakes and floodplains. This exchange is through the porous sediment, and it is characterized by the circulation of surface water into the alluvium and back to the river bed [6].

Moreover, fine sediments can move in suspension when the turbulent eddies have upward velocity components exceeding the fall velocity of fine sediments [7]. An increase in grain roughness can generate an increase in the vertical intensity of turbulence [8]. Moreover, strong upward turbulent ejections could provide the vertical anisotropy needed for suspension transport and the entrainment of the fine sediment [8–10]. In addition, the turbulent structures in smooth flumes have been studied by [11–14] and others, showing coherent structures such as individual hairpin vortices, which have a

small scale motion *l* ∼ *h*, hairpin vortices that have a large scale motion *l* ∼ 3*h*, and super streamwise vortices, which have a super scale motion *l* > 10*h*, where *l* is the streamwise scale of the vortices and *h* is the water depth. However, ejections in a rough-bed are smaller near the bed and the secondary currents can change the ejection patterns [15]. Furthermore, several researchers such as [12,16–19] have investigated turbulence characteristics considering a Cartesian plane with streamwise velocity fluctuations, *u* , and vertical velocity fluctuations, *w* , where the interaction in the second and fourth quadrant (ejections *u* < 0 & *w* > 0 and sweeps *u* > 0 & *w* < 0), respectively, are more frequent than the interactions in the first and third quadrants (outward interaction *u* > 0 & *w* > 0 and inward interaction *u* < 0 & *w* < 0), respectively, because of the mean shear stress is positive (i.e., *u w* < 0). Experimentally, researchers have used acoustic Doppler velocimetry (ADV), particle image velocimetry (PIV), laser Doppler velocimetry (LDA), or ultrasonic Doppler velocimeter (UDV) for the acquisition of velocity data. Niño and Musalem [16] used ADV for characterizing the turbulence interactions in a sand bed, reporting ejections and sweeps as the most frequent turbulence interactions. Niño and Musalem [16] also reported that the sweeps are more efficient for the entrainment of the particles into suspension than ejections. Sambrook Smith and Nicholas [17], Cooper et al. [18], and Chen et al. [20] implemented PIV to characterize the flow field and the turbulence properties. Sambrook Smith and Nicholas [17] experimented with the deposition of sand on gravel beds. They reported high velocities and high shear stresses that occur at the level of the crest of the major roughness elements and in their lee side; in addition, interactions such as sweeps and ejections are less frequent when roughness decreases, i.e., the main effect of fine sediment deposited in a gravel bed is to reduce the vertical velocities' gradients and shear stresses near the bed over the sand on the gravel bed.

The roughness and the flow depth can modify the turbulent structures at shallow flow conditions [21]. Roussinova [21] compared results for rough and smooth walls and found that the magnitude of the turbulence quantities are higher in the case of the rough wall, ejection events are prevalent over sweeps and in smooth wall cases, and there are ejections and sweeps in the vicinity of the hairpin vortex.

Manes et al. [22] compared the turbulence structure for permeable and impermeable beds in open channels, finding that large scale eddies generated within the surface flow have influence in the subsurface flow, and they think that it must be associated with pressure fluctuations.

Fourier series are often used to discuss the properties of a turbulent flow field [23,24]. The frequency analysis, for example, is derived from the Fourier spectrum. However, this analysis is for stationary signals and the Fourier transformation has no localization property, i.e., if a signal changes at one position, then the transform changes without the position of the change could be recognized "at a glance" [25,26]. For unsteady signals that have finite duration, such as in geophysical processes and hydrology, the wavelet transform and the cross wavelet transform are excellent tools for analyzing the physical relationships between the time series [27,28]. However, the open channel turbulent velocity fluctuations have been analyzed by Chen et al. [20] considering the wavelet coherency, i.e., measuring the wavelet correlation between two velocity series at a frequency *f* on a scale from zero to unity, finding that the wavelet analysis identified the scale of motions and the time of its occurrence.

The present study characterizes the turbulence structures associated with the interchange between subsurface and surface flows due to fine sediment, pumicite, deposited into the interstitial space, the pores of a gravel bed. Particle image velocimetry (PIV) was employed to measure the velocity field before and after the fine sediment was deposited. The velocity data were analyzed by a scatter plot of turbulent fluctuation, Fourier, and wavelet analysis, finding changes in the distribution of turbulence interactions for a flume with and without fine sediment deposits.

#### **2. Materials and Methods**

Experiments were carried out in an open channel with a sediment bed. With this experiment setup, it was possible to measure independently the surface and subsurface flow. In this research, an immobile

bed was considered. The fine material with which the bed is polluted is pumicite, or natural pozzolan, a raw material of minute grains of volcanic glass and ashes with the characteristics of clay.

#### *2.1. Experimental Set Up*

The experimental facility is an open channel, 0.03 m wide, 0.58 m long, and 0.63 m deep. The facility is divided into three parts. Upstream are the surface and subsurface input flows and the location of the seeding particles for PIV. At the center of the structure are the open channel, the bed, the point where the pumicite mixture is spilled, and the PIV measurement area. Downstream are the surface and subsurface output discharges (Figure 1a). The sediment bed has two layers of sediment. The surface layer is of gravel, 20 mm thick and median diameter *Dg* = 10 mm, and the subsurface layer is of fine gravel, 390 mm thick, and mean diameter of *Ds* = 2.45 mm (Figure 1b). The grain size distributions of sediment are shown in Figure 2. The density of both materials, gravel and fine gravel, is 2.65 g/cm3, whereas the pumicite has a characteristic diameter (*D*50) of *Dc* = 0.12 mm (Figure 2) and a density of 1.7 g/cm3. The pumicite was fed through an acrylic cone in the free surface, of 1.0 cm of diameter, during 6 s. The net weight of pumicite was 284.2 g for the experiments. Such feeding can simulate a soil failure that falls into the river. Furthermore, the pumicite is poured with a concentration of 57% by weight in water.

**Figure 1.** (**a**) experimental scheme used in the study; (**b**) experimental flume. The sediment column is used to measure high percorlations of fine sediment (pumicite) and subsurface flow. The surface flow rate is *Qsur* and subsurface flow is *Qsub*.

The measurement of flow rates was through volumetric gauges for the surface, *Qsur*, and subsurface, *Qsub*, flows in the experiments. A Photron FASTCAM Mini UX50 camera (San Diego, CA, USA) was used for Particle Image Velocimetry. This camera takes up to 2500 fps. A Nikon D3200 camera (Tokyo, Japan) was used for percolation analysis of pumicite. The Malvern Master Sizer 2000 equipment (Malvern, UK) of the Laboratory of Sedimentology of Universidad de Chile was used for measurement of the grain size distribution of pumicite. Experiments with several flow rates *Qsur* and *Qsub* are presented in Bustamante and Niño [29,30] for percolation of a variety of fine sediments, among them with pumicite. The experiments in this article are two, for the flow rates *Qsur* = 0.088 L/s and *Qsub* = 0.008 L/s. One with uniform flow with a friction slope *Sf* = 0.017, and water depth *H* = 67 mm, without fine sediment (E0) and the other of the same hydrodynamics conditions but with the spill of the fine sediment (E1). The article focuses, in this way, on the interaction between subsurface and surface flows due to the deposition of fine sediment in the bed by measuring turbulence in the surface flow.

**Figure 2.** The grain size distribution of sediments used in the study.

#### *2.2. Velocities*

PIV was implemented to measure the flow field. The PIV particle tracers were of diameter 0.12 mm and density 1.7 g/cm3. The flow was illuminated with a lamp with power of 18 W. The particle tracers are poured dry in low quantity in the stilling tank. Furthermore, for these conditions, the PIV particles are a good tracer because these particles only can move in suspension. Image acquisition was with the Fastcam Mini UX50 camera. It was implemented taking 250 fps during 140 s, before and after the spill of pumicite mixture, in the center of the cross section.

The velocity data were analyzed considering three methodologies. The first methodology is the estimation of the mean shear stress profile with the fluctuations of velocities with the spatially averaged open channel flow methodology proposed by Nikora et al. [31]. It was analyzed defining the total shear stress such as:

$$
\pi\_{\text{tot}} = \mu(\partial \overline{u}) / \partial z - \rho \left< \overline{u'w'} \right> - \rho \left< \overline{u} \overline{w} \right>\,,\tag{1}
$$

i.e., the mean shear stress has three components which are: viscous (*μ*(*∂u*)/*∂z*), turbulent (−*ρ u w* ), and form-induced (−*ρ u*˜*w*˜), where *μ* is the kinematic viscosity, *u* is the mean velocity, *u* and *w* are the velocity fluctuations in *x* (streamwise) and *z* (vertical) components, respectively. *u*˜ and *w*˜ are the form-induced disturbance in the flow variables (where *u*˜ = *u* − *u* and *w*˜ = *w* − *w*), *x* is the streamwise coordinate and *z* is the vertical coordinate. In addition, the spatial average mean velocity, < *u* > (*z*), results from this methodology. In Equation (1), the overbar means temporal average and the angular brackets mean spatial averages.

The second methodology considers the scatter plots of velocity fluctuations *u* and *w* , which were analyzed considering that they could be represented by an ellipse, such as Equation (2), where *α*, *Ra* and *Rb* are the angle of rotation, major axis and minor axis, respectively [32]. Parameter *α* is obtained as the slope of a linear fit for the scatter plot between *u* and *w* :

$$\frac{(\mathbf{x}\cdot\cos(a)+y\cdot\sin(a))^2}{Ra^2} + \frac{(\mathbf{x}\cdot\sin(a)+y\cdot\cos(a))^2}{Rb^2} = 1\tag{2}$$

Additionally, Niño and Musalem [16] and Wallace [19] measured the intensity of the ejections and sweeps in a measurement at one point by the parameter *K* = *u w* / < *u w* >. In this research, we have considered the parameter as: *K* = −*K* = −*u w* / | *u w* |. The percentage of distribution of the *K* parameter in an area of 73 mm × 68 mm at one point *x* and *z* of the experiment E0 is presented in Figure 3, where *K* <sup>90</sup> and *K* <sup>10</sup> are associated with the 10th and 90th percentiles of *K* , respectively. After finding the pairs (*u* , *w* ) associated with *K* <sup>90</sup> and *K* <sup>10</sup>, the *Ra* and *Rb* can be defined. The methodology used for characterizing the ellipses is presented in Figure 4.

**Figure 3.** Percentage of distribution of parameter K.

**Figure 4.** Methodology to characterize the scatter plot of velocity fluctuations through the ellipse parameters.

The third methodology was Fourier and wavelet analysis. The wavelet analysis has been implemented by [20,26,28,33,34] and others, in order to analyze different physical phenomena. In this research, the Morlet function was used as the mother wavelet. Morlet wavelet mother is a powerful tool for data analysis of low-oscillation functions [26].

#### **3. Results**

Hydraulic parameters for the experiment setup are the Froude number, *Fr* = *U*/ *gH*, Reynolds number, *Re* = *UH*/*ν*, *B*/*H* is the aspect ratio between wide flume (*B*) and depth water (*H*), *H*/*Dg* is the ratio of depth water and gravel diameter, *Sf* is the friction slope (bed slope), and *U* = *Qsur*/(*BH*), as Table 1 shows.


**Table 1.** Hydraulic parameters.

In experiment E1, the fine sediment in the channel was poured locally, at the beginning of the flume (*x* = 0.05 m), as a hyperconcentrated mixture of fine sediment with water. The dynamics of this mixture considers suspension, deposition, and percolation. After deposition, the fine sediment changes the permeability of the gravel bed. Thus, the interaction between surface and subsurface flow generates ejections of water from the bed to the water column. These jets eject subsurface water with fine sediment of low density deposited in the interstices of the gravel bed. The length of the jets has been approximately 15 mm, as shown in Figure 5a.

These sediment ejections generate a low entrainment of fine sediment into suspension transport or bedload transport, i.e., the sediment ejected is deposited in the neighborhood of the ejection, forming craters as bedforms (Figure 5c). Then, the experiment E1 is to characterize the ejection of the subsurface flow, which is modified by a lower permeability associated with the fine sediment and emerges to the surface current.

Figure 5 is localized in the tank at 20 cm downstream of the center part of the facility as shown in Figure 1. The velocity data are taken at the center of the cross section and at the point where the sediment ejections appear. Furthermore, the velocity fluctuations and wavelet analysis are analyzed at 1.5 cm up from the of pumicite level. The velocity profiles were analyzed considering spatially averaged open-channel flow proposed by Nikora et al. [35].

**Figure 5.** (**a**) sediment ejection for *Qsur* = 0.088 L/s y *Qsub* = 0.08 L/s; (**b**) measurement points; (**c**) top view of the sediment ejections.

#### *3.1. Velocities and Shear Velocities*

The velocity data were measured with PIV and processed with PIVLab. The spatial average mean velocity profile and the spatial average mean shear rate were obtained according to the double-averaged methodology proposed by Nikora et al. [35]. The spatial average mean velocity profile made dimensionless with the shear velocity for the experiments E0 and E1 are shown in Figure 6. Table 2 shows the shear velocities in both experiments, *u*∗*E*<sup>0</sup> and *u*∗*<sup>E</sup>*1. Shear velocities were calculated as *<sup>u</sup>*∗*Ei* = *τ*0*Ei*/*ρ*, where *τ*0*Ei* is the experimental shear stress in the bed and *i* = 0, 1, for experiment E0 or E1, respectively.

The dimensionless velocity profile for E1 is much more intense than for E0. This is because the shear velocity *u*∗*E*<sup>0</sup> is 39% higher than the shear velocity *u*∗*<sup>E</sup>*1. However, the ratio of the depth average mean velocities in E1 and E0, *UE*1/*UE*0, is 0.90. Therefore, the surface flow is faster in E0 than in E1 because the sediments ejection has a dominant vertical velocity. However, in an area without sediment ejection, the streamwise velocity in E1 has to be greater than E0 because there is outflow of the subsurface flow to the surface flow.

**Table 2.** Velocities and shear velocities.


**Figure 6.** Velocity profiles double-averaged, for both experiments E0 and E1, before and after the spill of pumicite mixture, respectively.

In Figure 5b are the points that were analyzed for the turbulence interactions in one sediment ejection—that is, upstream of the center of ejection (P1), the center of ejection (P2) and downstream of the center of ejection (P3). Turbulence fluctuations were analyzed under three approaches: scatter plot of *u* and *w* , velocity field, and wavelet analysis.

#### *3.2. Velocity Field*

Velocity vectors in the streamwise and vertical plane obtained through PVI processing are presented in Figure 7, for three different times, 40.00 s, 40.63 s, and 41.26 s. Figure 7 also presents the contour plot of the vertical velocity. The red polygon in that figure is limited where the vertical velocities are higher than the streamwise velocities in the ejection. Additionally, in Figure 7c, the blue lines show a coherent structure external to the movement we are observing, which is also a coherent structure from upstream. The vertical upward movement from the bed to the water column is associated with turbulent interactions of the ejection type. However, the sediment ejections reported in this research are associated with jets with sediment, due to the interaction between surface and subsurface flow, and this is different from the turbulent interactions of the ejection type investigated for rigid wall, as reported by [9,12,18,36,37] and others, who have analyzed turbulence interactions near the wall and have identified two main interactions as ejections and sweeps.

**Figure 7.** Velocity field for dimensionless vertical component (W/u\*) and vectors of streamwise velocity and vertical velocity, measured experimentally for E1, at **a** t = 40 s, **b** t = 40.63 s, **c** t = 41.26 s. Measured field was Δ*x* = L = 2.9 cm y Δ*z* = 3.5 cm and flow depth H = 6.7 cm. Direction of the flow: from right to left.

In E1, the vertical velocities, *w*, made dimensionless with mean streamwise velocity, *U*, were analyzed before and after the spill of pumicite mixture as a function of dimensionless time *tU*/*H*. In this case, the vertical velocity time series at the water depth for the three positions into the sediment ejection, i.e., in each position, P1, P2, and P3, shown in Figure 5b, the velocities series were taken in the entire water column. However, for the case of E0, with no-spill of pumicite mixture, the vertical velocities are those associated only to P1 and P2. The vertical in P1 is a measure of the vertical velocity at a point on the gravel ridges (Figures 7 and 8a), while the vertical in P2 is a measure in the gravel pores (Figures 7 and 8b). The vertical velocities in *E*0, for *z*/*H* > 0.1 and in positions P1 and P2 follow Taylor's frozen turbulence hypothesis in the streamwise direction quite well (Figure 8a,b). Vertical velocities are higher for *z*/*H* > 0.5 − 0.6 than near the bed at P1 or P2 (Figure 8a,b). For position P2 near the bed, for *z*/*H* < 0.1, there is an area with high vertical velocity for 11 < *tU*/*H* < 61 (Figure 8b) that is not detected for position P1. That is, the irregularities presented by the bed of gravel, with the ridges and low points, make Taylor's hypothesis of frozen turbulence invalid near the bottom.

**Figure 8.** Vertical velocity made dimensionless with *U*, *w*/*U* in the experiment E0, with no-spill of pumicite mixture, for: (**a**) upstream of the center of the ejection, P1, (**b**) center of the ejection, P2.

Conversely, in E1, after the spill of pumicite mixture, it is found that the ejections generate streamwise changes in the vertical velocities. At P2, in the gravel pore, an ejection take place, with high positive vertical velocities, for *z*/*H* < 0.2 and 0 < *tU*/*H* < 80 (Figure 9b). At P3, negative vertical velocities are predominant for *z*/*H* < 0.2. This behavior is associated with a current toward the bed (Figure 9c). At P1, vertical velocities are as the experiment E0 (Figure 9a). Thus, basically the image in time is as shown in Figure 7 in an instant, that is, the ejection is a quasi-steady feature of the flow.

**Figure 9.** Vertical velocity made dimensionless with *U*, *w*/*U* in the experiment E1 after the spill of pumicite mixture for: (**a**) upstream of the center of the ejection, P1; (**b**) center of the ejection, P2 and (**c**) downstream from the center of the ejection, P3.

#### *3.3. Scatter Plot of Velocity Fluctuations*

The turbulence interactions associated with fine sediment ejections (E1) were compared with turbulence interaction without fine sediment (E0). In addition, the experiments measured with PIV in this article (P1, P2 and P3) were compared with Acoustic Doppler Velocimetry (ADV) measurements of turbulent interactions for an open channel with a gravel bed [38]. Experimental setup for the acquisition of ADV data of [38] was a flow rate of 14 L/s, a gravel bed of 45 mm particle diameter, a flow depth *H* of 0.1 m and *hm*/*H* = 0.85, where *hm* was the location of the ADV with respect to the bed. The distribution of the turbulence fluctuations (*u* , *w* ) of an open channel in a Cartesian plane can be represented with an ellipse of negative slope and major axis (*Ra*) in the direction of quadrants 2 and 4 (*Q*2 − *Q*4) and minor axis (*Rb*) in direction of quadrants 1 and 3 (*Q*1 − *Q*3), as shown in Figure 10a,b, where the quadrants 1 to 4 are against clockwise circuit of the Cartesian plane and *Q*1 was (*u* > 0, *w* < 0). According to [12,16,17,19], the main turbulent coherent structures are ejections and sweeps, *Q*2 and *Q*4, respectively. The turbulent interactions in the narrow flume of this article have the same distribution as the widest flume of [38], i.e., an ellipse with a negative slope and the main turbulent coherent structures are ejections and sweeps; however, the magnitude of fluctuations are greater for [38], due to the turbulence is also greater, as shown Figure 10a,b. However, the pattern of those distributions change with respect to the turbulence in an open channel due to the presence of sediment ejections. These sediment ejections tend to increase the vertical fluctuations and decrease the

streamwise fluctuations. These changes generate that outwards and inwards interactions can become events with a greater probability of occurrence than in a regular open channel, i.e., the ellipse has a positive slope and a relationship between *Ra* and *Rb* close to 1 (Figure 7c–e).

**Figure 10.** Scatter plot for: (**a**) ADV data [38], (**b**) PIV data in experiment E0, (**c**) PIV data in experiment E1 after the spill of pumicite mixture for P1, (**d**) PIV data in experiment E1 after the spill of pumicite mixture for P2 and (**e**) PIV data in experiment E1 after the spill of pumicite mixture for P3.

#### *3.4. Frequency Spectra and Wavelet Analysis*

Power spectrum frequency for vertical and streamwise velocities fluctuation are calculated both for E0 and E1 in points P1, P2, and P3 (Figure 11). Power spectrum was made dimensionless with *u*2 <sup>∗</sup>*H*. In both cases, it is possible to identify the production zone and the inertial sub-range in the power spectrum for both cases before and after the spill of pumicite mixture. The inertial sub-range was considered between frequencies 10 and 20–30 Hz, where the slope of −5/3 is representative. The frequency of more than 20–30 Hz is noise in the PIV velocities. In the case of E0, in the production zone and in the three points of measurement, the energy is higher for streamwise velocity fluctuations than vertical velocity fluctuations. Conversely, for E1 after the spill of pumicite mixture, the vertical velocity fluctuation has higher energy in the production zone than the streamwise velocity fluctuation because the sediment ejection in P2 has dominant vertical velocities (Figure 11e). Then, according to the dynamics of sediment ejection (E1), where the vertical fluctuation is dominant, it can be seen that, at P2, the energy in the production zone for this components is highest, at P1 the energy in the production zone is the same as the experiment E0, and at P3 the energy in the production zone for the vertical component is large with respect to the streamwise component and can be considered as a downwelling point.

**Figure 11.** Power Spectrum Density for streamwise velocity fluctuation and vertical velocity fluctuations (*z*/*H* = 0.1), (**a**–**c**) in experiment E0, (**d**–**f**) in experiment E1, after the spill of fine sediment mixture.

Since the sediment ejection is a quasi-steady flow, then the local wavelet spectrum in the experiment E1, after the spill of fine sediment mixture, was implemented to analyze the *u* and *w* velocity components in the three points of the sediment ejection, P1, P2 and P3 (Figures 12 and 13). The wavelet spectra |*Wu*| <sup>2</sup> and |*Ww*| 2, corresponding to the wavelet spectra for *u* and *w* velocity fluctuations, respectively, were made dimensionless with dimensionless with *Hu*∗, as shown in Figures 12 and 13. The spectra have *λ*/*H* in vertical axes and *tU*/*H* in horizontal axes, where *λ* is the wavelength calculated as *λ* = *U*/ *f* , *f* is frequency and *tU*/*H* is the dimensionless time. Wavelet analysis allows for seeing the variations of the power spectrum over time, i.e., with this analysis, a turbulent structure is determined and the time the structure is present in the measurement time. In the streamwise direction, upstream of sediment ejection, P1, the energy is concentrated for *λ* = 10*H*, for a long period of time, 10 < *tU*/*H* < 82 (Figure 12). This wavelength is associated with a frequency of 0.06 Hz. According to the power spectrum shown in Figure 11e, it is a structure corresponding to a large-scale motion, i.e., this large-scale is present in almost all the time of the measurement time. In points P2 and P3, there is concentration of energy at this wavelength (frequency), but it is less intense than at P1, whereas, for *λ* = 5*H* and *λ* = 4*H*, there is a concentration of energy at all the measurement points, P1 to P3. This wavelength is associated with frequencies of 0.13 Hz and 0.16 Hz, corresponding to large-scale motions (Figure 11d–f). However, this structure is present only at certain points in time; the period most frequent is *tU*/*H* ∼ 50 (Figure 12a–c).

The wavelet spectrum for vertical fluctuation shows energy concentrations in P2 and P3 higher than those of P1; the highest concentration of energy is that of P2. This wavelet spectrum shows that, at point P2, during the entire measurement time, the vertical component of sediment ejection dominates, as shown in the spectrum of Figure 9b. According to the wavelet spectrum, Figure 13b,

*λ* = 10*H* has a high energy concentration for 10 < *tU*/*H* < 82. That wavelength corresponds to the frequency 0.06 Hz and is associated with a large scale of motion for P2 and P3 (Figure 11b,c). In addition, for *λ* = 5*H*, 4*H*, and 3*H*, there is a concentration of energy at points P2 and P3. This wavelength is associated with frequencies of 0.13 Hz, 0.16 Hz, and 0.21 Hz, respectively, corresponding to large-scale motions (Figure 11d–f). However, these structures are present only at certain points in time, the period most frequent is *tU*/*H* ∼ 40 (Figure 12a–c).

**Figure 12.** Wavelet spectrum for streamwise velocity fluctuations, *u* , in the experiment E1, after the spill of fine sediment mixture (**a**) P1, (**b**) P2, and (**c**) P3.

**Figure 13.** Wavelet spectrum for vertical velocity fluctuations, *w* , in the experiment E1, after the spill of fine sediment mixture, (**a**) P1, (**b**) P2, and (**c**) P3.

#### **4. Discussion**

The coherent structure in turbulence flows have been researched by [16,18,19,36,39] in open channels with permeable and impermeable beds, considering a steady flow. One methodology has been to analyze the Reynolds shear stress in quadrants of the Cartesian plane, with *u* and *w* and their distribution is described as an ellipse of negative slope and a turbulent event has a high probability to be in the quadrant *Q*2 and *Q*4, while *Q*1 and *Q*3 have a low probability of occurrence. In this study, we can see the same distribution of Reynolds shear stress in the case of open channel without fine sediment deposited into the bed (E0). However, the interaction between surface and subsurface flow in E1, after the spill of pumicite mixture, causes the analysis to change. The sediment ejection generates a quasi-steady flow from the bed toward the water column, where the vertical velocity component is higher than the streamwise velocity component, i.e., the turbulent interactions in *Q*1 and *Q*3 for that structure has a higher probability than in other research.

Furthermore, for isotropic turbulence in smooth open channel flow, Taylor's frozen turbulence hypothesis has been validated by [13,20,40] and others. Taylor's frozen turbulence hypothesis is also validated in experiment, E0, for *z*/*H* > 0.1. However, Ref. [41] recognizes the limitations of applying the Taylor's hypothesis. They considered an uncertainty for *z*/*H* < 0.1 because, in this region, the mean velocity and local velocity can diverge. Then, in experiment E0, the velocity difference between P1 and P2 for *z*/*H* < 0.1 is associated with a divergence of gravel pore velocities, so that, near the permeable bed, the Taylor's frozen turbulence hypothesis is invalidated, and these results are according to [41].

In experiment E1, fine sediment ejections are quasi-steady flows from the bed. Turbulence patterns change and turbulence becomes anisotropic because sediment ejection has a dominant vertical velocity. In this case, Taylor's frozen turbulence hypothesis is also invalidated.

The power spectrum density in open channels has been reported by [42–44]. They have validated the −5/3 slope in the measured spectra and show a greater spectral density of the streamwise component than the spanwise and vertical components. However, in this study, for the experiment E0, in the production zone, there is a greater spectral density of the streamwise component than the vertical component, but in the inertial subrange both components have a spectral density of similar magnitude. Conversely, in experiment E1, since the sediment ejection has a dominant upward movement, the spectral density changes. Upstream of the sediment ejection, at P1, in the production zone, the spectral density for streamwise component is greater than the vertical fluctuation component and the spectrum densities in both components are greater than the spectral density in experiment E0, whereas, in the center of the ejection, at P2, in the production zone and in the inertial subrange, the spectral density for vertical fluctuation component is greater than the spectral density for streamwise components. Downstream of the ejection in experiment E1, P3 can be considered as a downwelling point. It is important to note that sediment ejection increases the spectral energy, both in the production zone and in the inertial subrange.

Coherent structures of turbulence in open channels have been characterized by their sizes, such as small-scale motion, large-scale motion, and very-large scale motion by [13,20]. Streamwise scales of 3*H* and 10*H* allow coherent structures in these experiments to be classified as large-scale, i.e., hairpin vortex, and very large-scale motion, i.e., super streamwise vortex. Ref. [20] implemented the wavelet analysis to detect high concentrations of energy over time and their scale of motion, finding high concentrations of energy in *λ* = 3*H* and *λ* = 10*H*, i.e., they reported hairpin vortices and super stream vortices; however, they did not show the influence cone, so the presence of the super stream vortices is not clear. The large-scale motion and the very large-scale motions find in this research are not the same coherent structures reported by [13,20] because the pattern of the sediment ejection has a dominant vertical component. In addition, the energy concentration changes with the measurement position inside the sediment ejection.

The interaction between ejections, sweeps, outward, and inward interactions generates turbulent structures such as horseshoe vortices, hairpin vortices, shedding vortices, etc. Those turbulent structures can move sediment away from where they occurred, i.e., they can generate erosion or scour [45]. Interactions such as ejections and sweeps are more common in open channels; however, in the experiments presented in this article, this is not entirely true. Outward and inward interactions are more frequent than in open channels. Furthermore, pumicite is a fine cohesive material, so mechanisms such as particle rolling, sliding, and saltation were not observed in the present experiments. Therefore, no erosion or scour was observed due to the sediment expulsion. However, there is a constant interaction between the surface and subsurface flow. This dynamic is relevant considering that the fine material can be, for example, mining materials and, therefore, a constant exchange of toxic material can be generated from the hyporheic zone to the surface flow.

The experimental scale is small compared with natural environments, but the dimensionless parameters in the experiments presented in this article on the basis of the flow and sediment transport phenomena in mountain streams are believed to be correct. The dimensionless parameter *H*/*Dg* = 6.7 of the experiments is representative of macroroughness flow and there is no bedload transport, that is, *<sup>τ</sup>*<sup>∗</sup> <sup>=</sup> 0.00006 <sup>&</sup>lt; *<sup>τ</sup>*∗*<sup>c</sup>* <sup>=</sup> 0.035, where *<sup>τ</sup>*<sup>∗</sup> is the dimensionless shear stress, *<sup>τ</sup>*<sup>∗</sup> <sup>=</sup> *<sup>u</sup>*<sup>2</sup> <sup>∗</sup>/(*gRDg*), and *<sup>τ</sup>*∗*<sup>c</sup>* is the critical shear stress for the incipient motion of the sediment, with *R* = (*ρ<sup>s</sup>* − *ρ*)/*ρ*, *ρ<sup>s</sup>* and *ρ* are the sediment and water density, respectively [38]. However, the dimensionless parameter that is essential for this article is *ue*/*ws*, where *ue* is the entrainment velocity and *ws* is the fall velocity of the fine sediment [10]. The dimensionless parameter *ue*/*ws* > 1 would make entrainment possible. The flow average vertical velocity < *w* > in P2 is 0.007 m/s (Figure 9) and *ws* of the pumicite is 0.005 m/s, so, if *ue* is equal to < *w* >, < *w* > /*ws* > 1 and there is entrainment of the pumicite to the surface flow. The value of *w*/*U* in P2 observed in Figure 9 is in the order of 0.15, which shows that a significant part of the subsurface flow appears in the surface flow.

#### **5. Conclusions**

Fine sediment, such as pumicite, between the pores of a gravel river bed, can reduce the permeability and the initial porosity of the bed, modifying the roughness and hydraulic parameters of the subsurface flow. Pumicite has a low density, generating changes in the interaction between the subsurface and surface flow. These interactions are mass, momentum, and energy exchange, so the decrease of permeability of the gravel layer can generate an increase of vertical velocity and the turbulence intensities in the surface layer. Additionally, low velocities in the streamwise direction and high vertical velocities can break the streamwise structures associated with secondary currents near the bed and sediment ejections can be seen.

The sediment ejections change the patterns of turbulent structures and the distribution of the turbulence interactions, which means that the flow does not have a typical rough wall open channel flow turbulence. Additionally, the sediment ejections increase the energy both in the production zone and inertial subrange. Within the ejection, the vertical velocity component has the highest increase of energy in the center of the ejection. The sediment ejections could vary with the granular size of the subsurface layer and the density of the fine material, i.e., the low-density fine material in the subsurface layer encourages the presence of sediment ejections from the bed.

As future work, we will continue to evaluate the turbulent structures associated with sediment ejections in the presence of surface and subsurface flows. Fine sediments, with higher densities will be considered, for example, mining materials, such as tailings or metal concentrates, to evaluate the effect of particle density on the dynamics of the ejection. To characterize both spanwise and the sediment ejection in 3D, we will implement the technique Stereo PIV.

**Author Contributions:** Conceptualization, N.Y. and B.-P.N.; methodology, B.-P.N. and N.Y.; formal analysis, N.Y. and B.-P.N.; investigation, B.-P.N. and N.Y.; resources, N.Y.; writing—original draft preparation, B.-P.N.; writing—review and editing, N.Y.; supervision, N.Y.; project administration, N.Y.; funding acquisition, N.Y. In general, both authors contributed to the general development of the document through periodic meetings and research sessions. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors of this paper thank the financing of the Department of Civil Engineering, University of Chile, the Fondecyt Project 1140767, support from CONICYT through Beca Doctorado Nacional No 21181620 and Advanced Mining Technologic Center (AMTC) and CONICYT Project AFB180004.

**Acknowledgments:** The authors acknowledge Felipe Galaz for collaborating in conducting experiments.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**


c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Wavelet Coherency Structure in Open Channel Flow**

#### **Kebing Chen 1, Yifan Zhang <sup>1</sup> and Qiang Zhong 1,2,\***


Received: 4 July 2019; Accepted: 6 August 2019; Published: 12 August 2019

**Abstract:** Many studies based on single-point measurement data have demonstrated the impressive ability of wavelet coherency analysis to catch the coherent structures in the wall-bounded flows; however, the question of how the events found by the wavelet coherency analysis relate to the features of the coherent structures remains open. Time series of velocity fields in *x*–*y* plane of the steady open channel flow was obtained from a time-resolved particle imaging velocimetry system. The local wavelet spectrum found shows that one of the main energetic scales in open channel flows is 3*h* motions. The wavelet coherent coefficients of *u* and *v* series from the same measurement points successfully detected the occurrence and the scale of these 3*h* motions, and the phase angle indicates their inside velocity structure is organized by the Q2 and Q4 events. The wavelet coherency analysis between different measurement points further reveals the incline feature of the 3*h* scale motions. All the features of this 3*h* motion found by the wavelet coherency analysis coincide well with the flow field that is created by the passing of hairpin packets.

**Keywords:** wavelet coherency; Taylor's frozen turbulence hypothesis; scale; hairpin vortex packet; open channel flow

#### **1. Introduction**

The essential difficulties of turbulence are the numerous scale motions and the complex interactions between them. Therefore, the two major tasks of turbulent data analysis can be summarized: The first one is scale decomposition, which decomposes the different scale components from the original turbulent signal; the other one is the correlation of scales, which analyzes the interactions between different scales.

Traditional turbulent data analysis is based on the Fourier transform. Fourier analysis decomposes the turbulent signal into the sum of infinite sine and cosine functions with different frequencies. It can reveal the frequency characteristic of the turbulent signal, but it cannot determine when those frequency components occurred. This is due to the fundamental difference between the presuppositions of Fourier analysis and intermittency of the turbulent signal. Specifically, each frequency in Fourier analysis stands for one scale vortices, and the passing of one vortex will cause one period sine wave at the measurement point. The continuous sine waves from positive to negative infinity on the time axis make Fourier analysis implicitly presuppose that each scale vortices occupy the entire timeline. However, the vortices in turbulence occur intermittently, and they present in the portion of both time and space axis and only effect the local flow field.

Wavelet analysis overcomes this defect by using the waves fluctuating in the limited time as the basis functions, and the behavior of the signal at infinity does not play any role. Thus, wavelet analysis is performed in both time and frequency domain, allowing the definition of local spectral properties and the ability to zoom in on local features of the signal. Since Farge [1] introduced the wavelet analysis into turbulence, wavelets have become pervasive in turbulent signal analysis. Most of the applications use the wavelet analysis only as a scale decomposition tool, and the research on the correlation of scales are still based on the traditional correlation coefficient or the Fourier coherency analysis.

Following the definition of coherency in Fourier analysis, Liu [2] firstly defined the wavelet coherency and applied the signal from ocean wind waves. The wavelet coherency analysis can reveal the similarity of two signals at different scales and instants, which provides a powerful tool to reveal the correlation between two signals of different types or from different locations which is cased by the local coherent structures. The fundamentals of the cross-wavelet and wavelet coherency analysis were systematically introduced by Torrence and Compo [3], and Torrence and Webster [4] proposed a broadly applicable smoothing function for the wavelet coherency, which made wavelet coherency eventually become a universal data analysis method. Grinsted et al. [5] developed a software package that allows users to perform the cross wavelet transform and wavelet coherency analysis, and they applied the methods to the Arctic Oscillation index and the Baltic maximum sea ice extent record. Camussi et al. [6] applied the cross wavelet transform and wavelet coherency analysis to wall pressure fluctuation signals from a microphone pair in the incompressible turbulent boundary layers. The time instants corresponding to a local wavelet coherency overcoming the trigger level were selected as the coherent events, and conditional statistical analysis was performed on these selected events to show the physical features of these highly coherent events. Based on the conditional statistics, the highly coherent event was speculated as the footprint of the near-wall-sweep-type motion which are known to be closely associated with the presence of streamwise vortices embedded within the turbulent flow and located in the near-wall region.

From the previous literature review, it can be seen that most of the applications of the wavelet coherency was done on single point measurement signals. These single point signals contain the information of the events that happen in turbulence, but they cannot reflect the whole picture of the events, and the question of how these highly coherent events revealed by the wavelet coherent analysis from single point measurement data related to the two or three dimensional coherent structures in the flow field remains open. Although we can easily get two- or even three-dimensional flow fields in the laboratory today, a large number of tasks in the actual application scenario are still completed by single-point measurement. Therefore, it is instructive for the single point data analysis to explore the three-dimensional coherent structures corresponding to events in them. Present study employs the time-resolved particle imaging velocimetry system (TR-PIV) to get the two dimensional flow field time series, which makes it possible to find out the flow structures from the 2D flow fields at the corresponding wavelet coherent events from the single point time series. The following part of this paper is organized as follows. Section 1 describes the principle of the wavelet coherency analysis briefly. Section 2 introduces the experiment and data. Section 3 presents the results of wavelet spectral and coherency analysis. Section 4 summarizes major findings.

#### **2. Wavelet Coherency**

We limit ourselves to a brief introduction of the wavelet coherency. The detailed procedures can be found in [3–5]. The continuous wavelet transform of a velocity signal *u*(*t*) is defined as the convolution of *u*(*t*) with a scaled and normalized wavelet:

$$\mathcal{W}\_{\boldsymbol{u}}(\boldsymbol{a},t\_0) = \bigcap\_{-\infty}^{+\infty} [\boldsymbol{u}(t') \frac{1}{\sqrt{a}} \psi^\*(\frac{t\_0 - t'}{a})] \boldsymbol{d}t' \tag{1}$$

where ψ is the wavelet function; superscript \* denotes the complex conjugate; *a* and *t*<sup>0</sup> are the scale and position parameter, respectively; *Wu*(*a*, *t*0) is the wavelet transform coefficient of *u* at scale *a* and instant *t*0. The wavelet function ψ in the wavelet transform must have zero mean and be localized in both time and frequency domain [1]. Following Liu [2], Torrence and Compo [3] and Chui [7], *Water* **2019**, *11*, 1664

the complex Morlet wavelet is used in this study since it provides a good balance between time and frequency localization:

$$
\psi(t) = \pi^{-1/4} e^{i\omega\_0 t} e^{-\frac{\vec{\xi}^2}{2}} \tag{2}
$$

where ω<sup>0</sup> is a wavelet center frequency, here taken to be 6 to balance between time and frequency localization. From Equation (1), it can be seen that the wavelet function is stretched and shifted in time by varying the scale and position parameter. Thus, different scale motions at every instant can be separate out from the original signal by the wavelet transform. In order to facilitate comparison, the pseudo Fourier frequency corresponding to the scale *a* is usually determined by picking up the energy peak in the Fourier spectrum of the wavelet function. For the Morlet wavelet with ω<sup>0</sup> = 6, the Fourier frequency *f* is almost equal to the scale, e.g., Fourier frequency *f* <sup>0</sup> is 0.971Hz when *a* = 1, and the corresponding frequency *f* at scale *a* is determined by:

$$f = \frac{f\_0}{a} \tag{3}$$

When the pseudo Fourier frequency *f* is determined, the wave number *k* and wavelength λ can be obtained based on Taylor's hypothesis:

$$\begin{cases} \begin{array}{c} k = \frac{2\pi f}{L} \\ \lambda = \frac{2\pi}{k} \end{array} \end{cases} \tag{4}$$

where *U* is the mean streamwise velocity at that point.

Similar to Fourier analysis, the cross wavelet coefficient can be obtained:

$$\mathcal{W}\_{\rm uv}(a, t) = \mathcal{W}\_{\rm u}(a, t)\mathcal{W}\_{\rm v}^\*(a, t) \tag{5}$$

There is a Parseval relation [8]:

$$\int\_{-\infty}^{+\infty} [u(t)v^\*(t)]dt = \frac{1}{c\_\psi} \int\_{-\infty}^{+\infty} \int\_0^{+\infty} \mathcal{W}\_u(a,t)\mathcal{W}\_v^\*(a,t) \frac{da}{a^2} dt = \frac{1}{c\_\psi} \int\_{-\infty}^{+\infty} \int\_0^{+\infty} \mathcal{W}\_{\text{inv}}(a,t) \frac{da}{a^2} dt\tag{6}$$

where *c*<sup>ψ</sup> is a factor depended by the wavelet function:

$$\begin{cases} \mathcal{L}\_{\psi} = \int\_0^{\infty} \left| \hat{\psi}(\omega) \right|^2 \frac{d\omega}{|\omega|} \\ \quad \hat{\psi}(\omega) = \int\_{-\infty}^{+\infty} \psi(t) e^{-2i\pi\omega t} dt \end{cases} \tag{7}$$

In the turbulence situation, the velocity *u*(*t*) and *v*(*t*) are both real signal. Thus, we have:

$$\begin{split} E\_{\text{total}} &= \int\_{-\infty}^{+\infty} u(t)^2 dt \quad = \frac{1}{c\_{\psi}} \int\_{-\infty}^{+\infty} \int \left| \mathcal{W}\_{\text{uu}}(a\_{\tau}t) \right|^2 \frac{da}{a^2} dt \\ &= \frac{1}{c\_{\psi}k\_0} \int\_{-\infty}^{+\infty} \int \left| \mathcal{W}\_{\text{uu}}(\frac{f\_0}{f}, t) \right|^2 df dt \end{split} \tag{8}$$

where *E*total can be considered as the turbulent kinetic energy. Equation (8) shows that the wavelet transform coefficients actually represent the energy carried by the corresponding scale motions at the given instant. Thus *Wuu*(*f*, *<sup>t</sup>*) <sup>2</sup> is called local wavelet spectrum. The frequency spectrum can be defined as:

$$E\_{\rm un}(f) = \frac{1}{\mathcal{C}} \int\_{-\infty}^{+\infty} \left| \mathcal{W}\_{\rm un}(f, t) \right|^2 dt \tag{9}$$

*Water* **2019**, *11*, 1664

where *C* is the normalization constant which makes

$$\int\_0^{+\infty} E\_{uu}(f) df = u'^2 \tag{10}$$

in which *u*is the standard deviation of time series *u*(*t*).

For the correlation coefficient, we have the definition:

$$R\_{uv} = \frac{\int\_{-\infty}^{+\infty} [u(t)v(t)]dt}{\sqrt{\int\_{-\infty}^{+\infty} u(t)^2 dt} \cdot \sqrt{\int\_{-\infty}^{+\infty} v(t)^2 dt}} \tag{11}$$

Based on Equation (6), we can get:

$$R\_{uv}^2 = \frac{\Re\left[\int\_{-\infty}^{+\infty} \int\_0^\infty \mathcal{W}\_{uv}(a,t) \frac{dq}{a^2} dt\right]^2}{\left[\int\_{-\infty}^{+\infty} \int\_0^\cdot \left|\mathcal{W}\_u(a,t)\right|^2 \frac{dq}{a^2} dt\right] \cdot \left[\int\_{-\infty}^{+\infty} \int\_0^\cdot \left|\mathcal{W}\_v(a,t)\right|^2 \frac{dq}{a^2} dt\right]} \tag{12}$$

in which [·] is the real part of complex variables. Similar to the Equation (12), we can define the wavelet coherency as

$$\mathcal{R}\_{uv}^{2}(a,t) = \frac{\Re\left[\mathcal{W}\_{\text{uv}}(a,t)\right]^{2}}{\left|\mathcal{W}\_{\text{u}}(a,t)\right|^{2} \cdot \left|\mathcal{W}\_{\text{v}}(a,t)\right|^{2}}\tag{13}$$

However, the correlation between two signals at one instant obviously makes no sense. As noted by Liu [2], this coherency is identically one at all times and scales. This problem is circumvented by smoothing the wavelet coefficient along time or both time and scale axis before normalizing. The time and scale smoothing operator given by Torrence and Webster [4] is used in this study:

$$\begin{cases} \begin{aligned} \mathcal{R}\_{uv}^{2}(a,t) &= \frac{\left[\mathcal{S}\left[a^{-1}\mathcal{W}\_{uv}(a,t)\right]\right]^{2}}{\mathcal{S}\left[a^{-1}\left|\mathcal{W}\_{u}(a,t)\right|^{2}\right] \cdot \mathcal{S}\left[a^{-1}\left|\mathcal{W}\_{v}(a,t)\right|^{2}\right]} \\ \mathcal{S}\left[\mathcal{W}\right] &= \mathcal{S}\_{\text{scale}}(\mathcal{S}\_{\text{time}}(\mathcal{W})) \end{aligned} \tag{14} \\ \begin{aligned} \mathcal{S}\_{\text{scale}}(\mathcal{W}(a,t)) &= \mathcal{W}(a,t) \ast c\_{1} \Pi(0.6a) \\ \mathcal{S}\_{tim}(\mathcal{W}(a,t)) &= \mathcal{W}(a,t) \ast c\_{2} \varepsilon^{\frac{-\sigma^{2}}{2a^{2}}} \end{aligned} \tag{15} \end{cases} \tag{14}$$

in which the symbol \* denotes the convolution product; Π(0.6*a*) is a boxcar filter of width 0.6; *e* −*t* 2 2*a*2 is the absolute value of the Morlet wavelet; *c*<sup>1</sup> and *c*<sup>2</sup> are normalization coefficients to have a total weight of unity. The factor of 0.6 is the empirically determined scale decorrelation length for the Morlet wavelet. The wavelet-coherence phase difference is given by:

$$\phi(a,t) = \arctan\left\{\frac{\mathfrak{N}\left[\mathcal{S}\left[a^{-1}\mathcal{W}\_{\text{uv}}(a,t)\right]\right]}{\mathfrak{R}\left[\mathcal{S}\left[a^{-1}\mathcal{W}\_{\text{uv}}(a,t)\right]\right]}\right\} \tag{15}$$

in which [·] is the imaginary part of the complex variables.

The MATLAB function "wcoherence" is used to compute the wavelet coherency between two signals in this study. One can find the wavelet coherency of analytical signals in the documentation of the MATLAB function, which reveals the ability of wavelet coherency to analyze the scale and phase relations between two signals.

#### **3. Experiment**

Experiments were conducted in the Tsinghua tilting hydraulic flume, which is a closed-circuit open channel 20 m long and 0.3 m wide. The measurement section was set up 12 m downstream of the flume entrance. Eight ultrasonic water level sensors were set on the flume to monitor the water depth across the entire flume. The streamwise and wall-normal directions are denoted by *x* and *y*, respectively, and the corresponding components of fluctuating velocities are *u* and *v*. Velocity field measurements of the *x–y* plane in the middle of channel were made at a uniform flow condition as listed in Table 1. The water depth *h* is 2.9 cm. Thus, the flow in the central region can be considered as statistically two-dimensional [9,10]. The Reynolds number based on the section average velocity *U*m and water depth was 15,895 (*U*m was based on the discharge from the electromagnetic flowmeter), and the friction Reynolds number was 880 (*u*<sup>∗</sup> <sup>=</sup> *ghJ*).


**Table 1.** Flow condition and particle imaging velocimetry (PIV) parameters.

*J* = bed slope, *h* = water depth, ν = kinematic viscosity, *B* = channel width, *B*/*h* = aspect ratio, *Um* = the depth-averaged velocity, *u*\* = friction velocity, *Fr* = Froude number, *Re* = Reynolds number, *Re*τ = friction Reynolds number.

Instantaneous, two-dimensional velocity fields were measured in the streamwise-wall-normal plane (*x–y* plane) with a TR-PIV. The laser sheet was projected from the channel bed and it was located at the central line of the channel. The camera was set at the side of the channel and the laser sheet plane and the CMOS plane of the camera kept parallel to the mid-vertical plane of the channel. The PIV parameters are also listed in the Table 1. The exposure time was fixed at 150 μs as a compromise between minimizing image streaking and maximizing image lightness. The sampling frequency was 2500 Hz to obtain time-resolved series of 2D flow fields. A total number of 5596 images were obtained thus the time series contains 5595 continuous velocity fields. A greater sample number will be better but the capacity of the high-speed memory in the high-speed camera limits the number of images. Particle images were analyzed by using the iterative multi-grid image deformation method. Various test results of the PIV algorithm used in this paper can be found in the report of the 4th International PIV Challenge (the symbol of the PIV algorithm is TsU) [11]. The window size in the final iterative step is 16 × 16 pixels with a 50% overlap. A detailed description of the experiment system can be found in [12,13].

#### **4. Results**

#### *4.1. Wavelet Spectrum*

Figure 1a presents the wavenumber spectrum *Euu* for streamwise velocity fluctuation at point (*x*/*h* = 0, *y*/*h* = 0.5) from both Fourier and wavelet analysis. The biggest wavelength was chosen as 10*h* because the length of the time series is *tU*/*h* = 45.5, and this length is not enough for obtaining a creditable result for the larger scale motions. It can be seen from Figure 1a that the wavelet spectrum is much smoother and follow the Fourier spectrum well. This can be attributed to the fact that the wavelet spectrum presents an average of the Fourier spectrum weighted by the square of the Fourier transform of the analyzing wavelet shifted at wave number *k*, and it keeps the same power-law as in the Fourier spectrum [14]. Owing to the relatively small sample size, even the wavelet spectrum is still spiky. The TR-PIV system can capture the 2D flow field time series, but it is very hard to get big sample sizes because of the storage and transport difficulties of the huge amount of data. However, the main goal of this study is to show the physical meaning revealed by the wavelet analysis on the turbulent signal from open channel flows instead of presenting accurate wavenumber spectra.

Figure 1b shows the pre-multiplied power spectrum, *kEuu(k)*, at point (*x*/*h* = 0, *y*/*h* = 0.5) from the wavelet spectrum. Based on Equation (10), when the horizontal axis is set to logarithmic coordinates, we have

**Figure 1.** The wavenumber spectrum (**a**), pre-multiplied spectrum (**b**) and local wavelet spectrum (**c**) of the streamwise velocity series at the measurement point (*x*/*h* = 0, *y*/*h* = 0.5). The black solid straight line in (**a**) shows the classical −5/3 law in the turbulence. The black dash line in (**b**) and red dash line in (**c**) is marked the strongest scale in the spectrum.

Therefore, the whole area under the pre-multiplied spectrum curve in the semi-log plot is directly related to the value of turbulence intensity, and the area under a small section of the spectrum curve can be considered as the strength of the corresponding scale motions. The strongest scale in Figure 1b is approximately 3*h* as marked by the black dash line, which means the 3*h* scale motions are one of the main energetic structures in the outer layer. In the scales greater than 3*h*, there may exist another energetic mode around 10*h*. This two-energetic-mode feature of the pre-multiplied spectrum is similar to the results from other wall-bounded flows [15–17]. By the coherent structure classification in open channel flow [12,13], the *h* and 10*h* order motions can be classified as large- and very-large-scale structures.

The traditional spectrum gives the general information about the energetic scales but cannot show the time instants when the local event happens. The contour map of the local wavelet spectrum, |*Wu*| <sup>2</sup> is shown in Figure 1c. There are two high energy regions at scale approximately 3*h* and 10*h*, respectively, which coincides with the result in pre-multiplied spectrum in Figure 1b. From |*Wu*| 2, we can see that the strongest 3*h* and 10*h* scale motions pass the measurement point during 17 < *tU*/*h* < 21 and 25 < *tU*/*h* < 35, respectively.

In order to show the entire flow structures when the local wavelet spectrum shows high energy, fluctuating velocity field series from *tU*/*h* = 11 to 35 are pieced together based on Taylor's frozen turbulence hypothesis in Figure 2, following Zhong et al. [12].

**Figure 2.** Fluctuating velocity fields pieced together based on Taylor's frozen turbulence hypothesis. The red arrows indicate the strong Q2 and Q4 events in the flow field, and the red solid lines indicate the inclined shear layers.

The convection velocity in Taylor's frozen turbulence hypothesis is the mean velocity of the whole field. The major characteristic of flow field during 17 < *tU*/*h* < 21 is an inclined shear layer (marked by the red solid lines in Figure 2) with the strong ejections (marked by the red solid arrows in Figure 2). The ejections are also called the Q2 event, where the fluctuating streamwise velocity *u* < 0 and vertical velocity *v* > 0. When *u* > 0 and vertical velocity *v* < 0, there is a sweep motion, or Q4 event. The streamwise scale of this inclined shear layer structure is approximately 2.5*h*. The features of this inclined shear layer structure are similar to the flow field in the reference [12,18–20], and it is usually considered as the typical sign of the passing of a hairpin packet [19,21]. The transport effect of the hairpin vortices in the packet induce the strong Q2 events, and the shear line is inclined because the heads of hairpin vortices in the packet describe an envelope inclined at 15–20◦ with respect to the wall. Considering the information from flow field and the results in previous literature, the highly energetic region at 3*h* scale in the wavelet spectrum reveals the passing of the hairpin packet. In Figure 2 it can be seen that from *tU*/*h* = 11 to 30 there are at least six typical inclined shear layers located almost end to end, and the wavelet spectrum indeed shows high energy in the corresponding duration and scale.

For the 10*h* scale highly energetic region, the flow field during 25 < *tU*/*h* < 35 shows hairpin packets firstly and then strong Q4 events (marked by the red dash arrows in Figure 2). It indicates that the 10*h* order motions consist of two different types of smaller structures, hairpin packets, and Q4 events. This agrees with the phenomena reported by Zhong et al. [12], Adrian and Marusic [22], Zhong et al. [12] and Zhong et al. [13] suggested a super-streamwise vortex model for the 10*h* order motions in open channel flows. The super-streamwise vortices rotate around the *x* axis. The strong Q2 events from hairpin packets constitute the upward movement and the Q4 events are the downward movement of the super-streamwise vortices.

From the above discussion, the energetic scales can be presented by the traditional spectrum analysis, and the wavelet spectrum can show not only the energetic scales but also the moment these energetic structures occur. However, wavelet spectrum contains no information about the organization of these structures. The following analysis will show that the wavelet coherency can reveal more details about the energetic structures.

#### *4.2. Wavelet Coherency*

The wavelet coherency between time series *u* and *v* at the same measurement point (*x*/*h* = 0, *y*/*h* = 0.5) is shown in Figure 3. The advection velocity *U* is the mean velocity of the whole field. The white dash line marks the cone of influence of the boundaries. The biggest wavelength was chosen as 10*h* as in Figure 1. It can be seen from Figure 3 that the most remarkably coherence appears in 13 < *tU*/*h* < 32 (marked by the white rectangle) near the scale 3*h* (marked by the red dash line) and in *tU*/*h* >35 at the scale 3*h* to 10*h*. Most of the coherence area of *tU*/*h* > 35 is beyond the white dash line, which is the cone of influence for the wavelet, and the wavelet coherent value is untrusted because of too close to the beginning or ending instantaneous. Thus, we only focus on the 13 < *tU*/*h* < 32 region. The phase angle in this area is about π. This highly coherent area indicates that there are 3*h* scale structures continuously passing the measurement point during 13 < *tU*/*h* < 32. From Figure 2, one can see that there are several hairpin packets in the flow field during 13 < *tU*/*h* < 32. Thus, the highly coherent area near the scale λ/*h* = 3 in Figure 3 reveals the passing of several hairpin packets. In addition, the π phase angle further reveals the velocity structure inside the hairpin packets. As shown in Figure 4, when the hairpin packet passes the measurement point, the Q2 events lead to that the streamwise fluctuation *u* appearing as positive firstly, and the vertical fluctuation *v* appearing as negative. Both *u* and *v* cross 0 when the shear layer passes the measurement point, and then *u* and *v* turn into negative and positive, respectively. Therefore, the phase difference between time series *u* and *v* has to be approximately 180◦ when hairpin packets are passing. The π phase angle reveals the feature of the hairpin packets that the Q2 and Q4 events dominate their inside velocity structures.

**Figure 3.** Wavelet coherent coefficients between the streamwise and vertical velocity series at the measurement point (*x*/*h* = 0, *y*/*h* = 0.5). The white box marks the time duration as the same as that in Figure 2, and the red dash line indicates the 3*h* scale. The high coherent area appears in the white box and the phase difference is approximately π.

**Figure 4.** The explanation of the phase angle in the wavelet coherency between the streamwise and vertical velocity series at the same measurement point when the incline structures pass. When the hairpin packet passes the measurement point, the Q2 events lead to that the streamwise fluctuation *u* appears in positive firstly, and the vertical fluctuation *v* appears in negative. After the inclined shear layer passes the measurement point, *u* and *v* turn into negative and positive, respectively.

The *u* series at (*x*/*h* = 0, *y*/*h* = 0.5) was chosen as the fixed point to do the wavelet coherency with other three measurement points. Figure 5 shows the wavelet coherency of *u* series at (*x*/*h* = −0.35, *y*/*h* = 0.15), (*x*/*h* = 0, *y*/*h* = 0.15) and (*x*/*h* = 0.35, *y*/*h* = 0.15), respectively. Points (*x*/*h* = −0.35, *y*/*h* = 0.15) and (*x*/*h* = 0.35, *y*/*h* = 0.15) are located at the upstream and downstream of point (*x*/*h* = 0, *y*/*h* = 0.5), respectively, and (*x*/*h* = 0, *y*/*h* = 0.15) is directly below the fixed point. From Figure 5a, there are two highly coherent areas during 11 < *tU*/*h* = 35 near the scale 3*h*. The first one is roughly from *tU*/*h* = 11 to 23 (marked by the white box). The second one is approximately from *tU*/*h* = 27 to 30 and much weaker than the correlation during the same duration in Figure 3, which can be attributed to the fact that the correlation reduces when the distance between two measurement points increases. Referring to Figure 2, these two highly coherent areas are related to the typical hairpin packets during 11 < *tU*/*h* < 21 and 27 < *tU*/*h* < 29, respectively.

**Figure 5.** Wavelet coherent coefficients between the streamwise velocity series at different measurement points. The fixed point is (*x*/*h* = 0, *y*/*h* = 0.5), the moving points are (*x*/*h* = −0.35, *y*/*h* = 0.15), (*x*/*h* = 0, *y*/*h* = 0.15), and (*x*/*h* = 0.35, *y*/*h* = 0.15). The highly coherent area reduces and the phase angle increases when the lower point moving from the upstream to downstream of the fixed upper point.

The comparison of the coherent area during 11 < *tU*/*h* < 21 between three points in Figure 5a–c shows the coherent area reduces when the measurement point moves from the upstream to downstream. It means the overlap time duration of the 3*h* scale motions between fixed and moving points reduces. This is caused by the incline feature of the 3*h* scale motions. Figure 6 sketches cartoons to show the explanation. Two solid red dots represent two measurement points. Blue and yellow shapes stand for the inclined 3*h* scale structures. The blue and yellow shape mark the starting and end time instant, respectively, that both measurement points are inside the structure. *D* is the distance between the centers of the blue and yellow shapes. Thus *D*/*U* is the time interval of both measurement points located in the structure, which is the highly coherent area in Figure 5. It can be seen from Figure 6 that *D* is the greatest when the lower point is located at the upstream of the upper point (corresponding to

Figure 5a). *D* decreases when the lower point moves from the upstream to downstream, as shown in Figure 5b,c, because of the incline feature of the hairpin packets. In addition, the phase angle increases when the lower point moves from the upstream to downstream. This is also cased by the incline feature of the hairpin packets. When the lower point is located at the proper location of the upstream of the upper point, the two points almost enter the structure at the same time, as shown in Figure 6a. Thus, the typical signal of the hairpin packet presents at the same time in the velocity series in the two points and the phase difference is small. When the two points are on the vertical line, as shown in Figure 6b, the lower point enters the structure after the upper one enters the structure for a while. Therefore, there is hysteresis between the phases of the signal of the corresponding scale. As the lower point keeps on moving downstream, the phase difference will increase.

From the above discussion, the wavelet coherency analysis between the streamwise and wall-normal velocity series at the same point and the streamwise velocity series from different points not only detects the occurrence and scales of hairpin packets in open channel flows, but also reveal the internal velocity organizations and the tilt geometry of hairpin packets. The wavelet coherent coefficient and the phase angles present a powerful tool for the detection of the energetic coherent structures and the analysis of their internal organizations.

**Figure 6.** The explanation of the reducing wavelet coherent area and the increasing phase angle when the lower point moving from the upstream to downstream of the fixed upper point. When the lower point is located at the proper location of the upstream of the upper point, the two points almost enter the structure at the same time, as shown in (**a**). Thus, the typical signal of the hairpin packet presents at the same time in the velocity series in the two points and the phase difference is small, as shown in Figure 5a. When the two points are on the vertical line as shown in (**b**), the lower point enters the structure after the upper one enters the structure for a while. There is hysteresis between the phases of the signal of the corresponding scale (see Figure 5b). As the lower point keeps on moving downstream, the phase difference will increase, as in (**c**).

#### **5. Concluding Remarks**

Many studies based on single-point measurement data have demonstrated the impressive ability of the wavelet coherency analysis to catch the coherent structures in the wall-bounded flows, however, the question that how the events found by the wavelet coherency analysis based on single point measurement data relate to the features of the three-dimensional coherent structures remains open. Present study employed the TR-PIV system to produce a time series of a two dimensional flow field in the streamwise-wall-normal plane of steady open channel flows, which makes it possible to find out entire flow structures from the 2D flow field of the corresponding wavelet coherent events. Wavelet coherency analysis was applied on the velocity series from a single or multi-point to find the highly energetic and coherent events. The fluctuating velocity fields during the time period with high energy and high wavelet coherence were pieced together based on Taylor's frozen turbulence hypothesis to

explore the corresponding entire coherent structures of the wavelet events. The major finding of this study are as follows:


**Author Contributions:** Conceptualization, Q.Z.; methodology, K.C. and Y.Z; data curation, Y.Z.; writing—original draft preparation, K.C.; Y.Z., writing—review and editing, Q.Z.; supervision, Q.Z.; funding acquisition, Q.Z.

**Funding:** The study was supported by the National Natural Science Foundation of China (Grant No.51809268) and the Fundamental Research Funds for the Central Universities (China Agricultural University, Project No. 2019TC043).

**Conflicts of Interest:** The authors declare no conflict of interest.


#### **Symbol Table**

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Bedform Morphology in the Area of the Confluence of the Negro and Solimões-Amazon Rivers, Brazil**

#### **Carlo Gualtieri 1,\*, Ivo Martone 1, Naziano Pantoja Filizola Junior <sup>2</sup> and Marco Ianniruberto <sup>3</sup>**


Received: 6 May 2020; Accepted: 3 June 2020; Published: 6 June 2020

**Abstract:** Confluences are common components of all riverine systems, characterized by converging flow streamlines and the mixing of separate flows. The fluid dynamics of confluences possesses a highly complex structure with several common types of flow features observed. A field study was recently conducted in the area of the confluence of the Negro and Solimões/Amazon Rivers, Brazil, collecting a series of Acoustic Doppler Current Profiler (ADCP) transects in different flow conditions. These data were used to investigate the morphology of the bedforms observed in that area. First, the bedforms were mostly classified as large and very large dunes according to Ashley et al. (1990), with an observed maximum wavelength and wave height of 350 and 12 m, respectively. Second, a comparison between low flow and relatively high flow conditions showed that wavelength and wave height increased as the river discharge increased in agreement with previous literature studies. Third, the lee side angle was consistently below 10◦, with an average value of about 3.0◦, without flow separation confirming past findings on low-angle dunes. Finally, a comparison between the bedform sizes and past literature studies on large rivers suggested that while several dunes were in equilibrium with the flow, several largest bedforms were found to be probably adapting to discharge changes in the river.

**Keywords:** river hydrodynamics; ADCP; bedforms morphology; river confluence; Amazon River

#### **1. Introduction**

Bedforms are a very common feature in alluvial fluvial channels as a result of the unstable interaction between water flow, sediment transport and bed morphology [1,2]. Bedform analyses range from fundamental analytical descriptions [3], to detailed numerical simulations [4–6] and laboratory measurements [7–11], to large-scale field measurements [12–15]. In the engineering context, these analyses are fundamental in predicting discharge, flow resistance [16] and sediment transport in rivers. In geology, bedforms are studied as primary sedimentary structures that are forming at the time of deposition of the sediment and reflect the characteristics of the depositional environment [17,18].

Bedforms are often approximated by triangular shapes, but natural bedforms present a more complicated bi and three-dimensional morphology. The broadest classification of unidirectional flow bedforms is based on the flow regime under which the bedforms develop. Simons and Richardson [19] distinguished a lower flow regime, for subcritical flows, associated with ripples and dunes, and through a transitional flat bed, an upper flow regime, for supercritical flows, with standing waves and antidunes. The upper flow regime is uncommon in deep rivers, where the observed bedforms are mainly ripples and dunes [13,20]. Several criteria were proposed to identify ripples and dunes. Ripples range in length from approximately 0.05 to about 0.6 m and in height from 0.005 to just less than 0.05 m [17],

but, according to Yalin [21], ripples have length and height less than 0.6 and 0.04 m, respectively. Ashley et al. [22] classified dunes into four groups based on height (*Hbf*) and wavelength (λ*bf*), as listed in Table 1.


**Table 1.** Classification of bedforms according to Ashley et al. [22].

While ripple size is scaling with the size of the grains of the bed and it is independent of the flow depth [23], dunes size is related to flow conditions. Several equations were proposed in the literature to predict dune sizes and their adjustment to flow changes [9–11,13,23–26]. Dune development is controlled by the interaction among the flow, sediment transport, and dune form. Deformation and adaptation have been recognised as the main mechanisms responsible for dune development [9]. Even where the statistical descriptors of the dune population have converged (equilibrium state) and the reach-averaged bed shear stress is constant, variability in dune shape is caused by continuously deformation as they migrate. Additional variability is introduced by the adaptation of dunes to changes in flow. Dune adaptation is ubiquitous because river flow is typically both unsteady and non-uniform at the temporal and spatial scales that are needed for dunes equilibrium [9]. However, equilibrium conditions are not frequent because they require time and a sufficient rate of sediment transport. Bedforms increasing/decreasing in size over the time are called developing/diminishing dunes, while they may be considered in equilibrium if they migrate without any change in shape or mass, i.e., deformation [20]. The temporal lag of the development of dunes relative to their formative flow, i.e., dune hysteresis [9], directly depends on their size [13,27]. Dunes are generally seen to grow in size during the rising stage, reaching their maximum development after peak discharge. However, it was found that water depth and flow velocity have separate effects on dune adaptation. Dunes crests/troughs do not respond simultaneously to changes in flow and crest flattening, i.e., increasing bedform length while height is decreasing, is related to decreasing depth and increasing flow velocity [9]. During the falling stage, dunes may stretch and flatten with a rapid downstream migration of the lee side, remaining as large bedforms during the subsequent low-flow period [9,13,20].

According to Best [1], the fluid dynamics of asymmetric river dunes with an angle-of-repose lee side and generated in a steady, uniform unidirectional flow, is characterised by (1) accelerating flow over the dune stoss side; (2) flow separation or deceleration in the lee of the dune, with reattachment at from 4 to 6 dune heights downstream; (3) a shear layer bounding the separation zone, which divides this recirculating flow from the free stream fluid above; (4) an expanding flow region in the dune lee side; and (5) downstream of the reattachment point, a new boundary layer that grows beneath the wake along the stoss slope of the next dune downstream.

Such flow structure has many important implications for flow resistance [1]. The differential pressures generated by flow separation and flow acceleration/deceleration associated with the dune form generate a net force on the dune, called the 'form drag', which, together with the grain roughness drag, called the "skin drag", determines dune morphology and flow resistance [28,29]. However, low-angle dunes (lee-side angle < 10◦) do not possess a zone of permanent flow separation, and those with lee-side angles < 4◦ are believed to possess no flow separation at all, resulting in lower energy losses [15]. In this regard, a recent analysis of high-resolution bathymetry data demonstrated that the largest rivers on Earth are characterized by low-angle lee-sides (mean ~10◦) [15].

Dune features are expected to show an even greater complexity at river confluences, as these are characterized by very complex hydrodynamics and morphodynamics located in the Confluence Hydrodynamic Zone (CHZ) [30,31]. The CHZ is often characterised as a stagnation zone, a velocity deflection and re-alignment zone, a separation region with recirculation, a maximum velocity and flow recovery region [32]. The central part of the CHZ ends where flow recovery starts, while the CHZ ends where the flow is no more significantly affected by the confluence. The hydrodynamics and morphodynamics within the CHZ are influenced by: the planform of the confluence; the momentum flux ratio of merging streams; the level of concordance between channel beds at the confluence entrance; and differences in the water characteristics (temperature, conductivity, suspended sediment concentration) between the incoming tributary flows lead to the development of a mixing interface and may impact local processes about the confluence [32,33]. Confluence bed morphology is characterized by the presence of a scour hole, bars (tributary-mouth, mid-channel and bank-attached bars), and a region of sediment accumulation near the upstream junction corner [30].

This paper presents and analyses the morphology of the bedforms observed during two field surveys carried out at the confluence between Rio Negro and Rio Solimões in the Amazon Basin. The surveys were carried out using Acoustic Doppler Current Profiler (ADCP) during low flow in 2014 and in relatively high flow in 2015. The objectives of this study are to (1) describe the bedforms characteristics at the Negro/Solimões confluence; (2) compare bedform characteristics in different flow conditions; and (3) compare bedform scales with those derived from literature theoretical/empirical equations and from past field studies conducted in large rivers.

#### **2. Field Site and Campaigns. Basic Results on Hydrodynamics and Morphodynamics**

#### *2.1. Hydrological and Sedimentological Background*

The study area is centered about the confluence of the Negro and Solimões Rivers, located near Manaus in Northern Brazil, where these rivers merge to form the Amazon River, (Figure 1). The Negro and Solimões confluence ranks among the largest on Earth and is famous for visibly revealing the meeting of the black (Negro) and white (Solimões) waters of the two rivers.

**Figure 1.** Map of area of the Negro/Solimões confluence. White and coloured solid lines represent transects surveyed with Acoustic Doppler Current Profiler (ADCP) during FS-CNS1 and FS-CNS2 (Long1: red line; N-CNS: blue line; S-CNS: green line). Background: Landsat 8 image. Map datum: WGS84; projection: UTM 21S.

The distinct waters of these two rivers are related to the geology (soils, soil coverage, etc.) of the two respective catchments within the Amazon Basin. The Negro sub-basin is located in the North draining the western slopes of the Guyana Shield, which is characterised by gentle gradients and densely vegetated margin, which, in turn, implies a low sediment production [34], while the Solimões catchment includes the eastern margin of the Andes, where the combination of high declivity and erodible rocks gives origin to high sediment production [35,36]. The mean water discharge of the Negro and Solimões Rivers is about 30,000 and 100,000 m3/s, accounting for 14% and 49%, respectively, of the total freshwater discharge of the Amazon River into the Atlantic Ocean [34]. The two rivers have a different hydrological cycle: the Negro has two distinct discharge peaks along the year, the first of low-amplitude during the first three months of the year, and the second larger in the middle of the year; the Solimões has one peak between May and June [34]. Hydrologic data collected at the fluviometric stations, located in Tatu-Paricatuba (Negro River) and Manacapuru (Solimões River) (Figure 2), (www.ore-hybam.org/) show that low-flow conditions are occurring during the local Autumn, while the high-flow conditions occur during the local Winter season (June–July) (Figure 3). In terms of sediment load, the difference between the two rivers is even larger. The Solimões River has at the Manacapuru station an average load of suspended solids of 14,174 Kg/s, accounting for more than half of the total load of the Amazon River into the Atlantic Ocean, while the Negro River has at Paricatuba station an average suspended load of 254 Kg/s [37].

**Figure 2.** Location of the Tatu-Paricatuba and Manacapuru stations.

**Figure 3.** Discharge time series from 2008 to 2016.

#### *2.2. Field Campaigns*

The field campaigns were carried out about the Negro/Solimões confluence in the Amazon River basin within the EU-funded Clim-Amazon Project, to study low flow (October 2014, FS−CNS1 campaign) and relatively high flow (April/May 2015, FS−CNS2 campaign) conditions [38–43]. The discharges measured during the campaigns were within the typical observed values for the seasons and daily differences were at most 5%. The surveys were carried out using an Acoustic Doppler Current Profiler (ADCP) as well as a multi-parameter probe for the measurement of water physico-chemical parameters (temperature, conductivity, turbidity, etc.) and total suspended sediment (TSS) concentration. During the surveys, a Teledyne RDI 600 kHz Rio Grande ADCP was used to collect flow velocity, water column backscattering [44] and water depth at key locations about the confluence. In total, 98 cross-sectional transects were collected. In addition, 2 and 3 longitudinal profiles along both sides of the Amazon River were collected in FS−CNS1 and FS−CNS2 surveys, respectively.

#### *2.3. Basic Observations of Hydrodynamics and Morphodynamics about the Confluence*

Table 2 lists the main flow properties of the Negro and Solimões rivers measured just upstream of the confluence (N-CNS and S-CNS, three for each river and each field campaign) during the two surveys. Large differences in discharge and flow velocities were observed in the Solimões River between the two surveys, whereas, on the Negro River, these differences were smaller.

**Table 2.** Main flow properties of Negro and Solimões Rivers during FS−CNS1/FS−CNS2/FS−CNS3. *Q* = median discharge; *A* = median cross-sectional area; *W* = channel median width; *hmed* = median depth; *W*/*hrect* = median of the aspect ratio; *Vavg* = median of the cross-section velocity (*Q*/*A*); *Vdepth-avg* = median of the depth-averaged velocity; *Dir* = median of flow direction degrees from North; *Vmax* = maximum depth-averaged velocity.


From FS−CNS1 to FS−CNS2, the maximum depth-averaged velocity was almost constant in the Negro River, but the Solimões River increased from 2.2 to 2.6 m/s. Furthermore, from low to high flow conditions, the Negro channel increased in depth, from 24 to 31 m, but not in width, whereas in the Solimões River the width increased from 1.6 to 1.9 km and the depth from 27 to 28 m. Finally, from FS−CNS1 to FS−CNS2, the median flow direction in the Negro River remained unchanged, whereas in the Solimões River a significant change in direction occurred. It is worth noting that the confluence junction angle is of about 65◦ (Figure 1). At the confluence entrance, the Negro channel is almost aligned with the Amazon channel, while the Solimões-Amazon waters must undergo a large change in flow direction (60◦–70◦) to enter in the Amazon channel. Common hydrodynamic features [32] already noted in past confluence studies were present even about the Negro/Solimões confluence. The approximate location of those features is shown in Figure 4, where they are numbered as: (1) the stagnation zone; (2) the region of deflection; (3) the region of maximum velocity; (4) the downstream separation zone; (5) the beginning of the region of flow recovery; and (6) the end of the CHZ. Further details about hydrodynamics are described in [41].

Ianniruberto et al. [43] identified that: upstream of the confluence, the Negro side is mostly characterized by a rocky bed with fine sand cover, whilst on the Solimões side the river bed consists predominantly of sand, and a sediment deposition region occurs at the junction corner in correspondence with the stagnation zone (Figure 4, feature 1); the central part of the CHZ is characterised by a sediment by-pass region exposing eroded Cretaceous bedrock, with a scour hole corresponding to the region of maximum velocity (Figure 4, feature 3); deposited sediment forming a bank-attached bar on the Solimões side of the CHZ (Figure 4, feature 4); avalanche faces of sediments at the Solimoes mouth; a bedrock terrace was observed towards the downstream end of the CHZ (Figure 4, feature 5), marking the transition from rocky to alluvial bed, where the bedforms were found.

**Figure 4.** Map of the depth-averaged velocities about Negro/Solimões confluence on (**a**) FS−CNS1 and (**b**) FS−CNS2, with the location of the hydrodynamics features. Legend: (1) stagnation zone; (2) region of deflection; (3) region of maximum velocity; (4) downstream separation zone; (5) beginning of the region of flow recovery; and (6) end of the confluence hydrodynamic zone.

#### **3. Bedforms morphology. Results and Discussion**

#### *3.1. Hydrodynamics and Sediment Transport Parameters along the Longitudinal Transects. Results*

The location and length of the two longitudinal transects collected in the area of the confluence during FS−CNS1 and FS−CNS2, approximately 900 m away from the Solimões/Amazon right bank, are presented in Table 3 and Figure 1.


**Table 3.** Locations and length of the ADCP transects collected in FS−CNS1 and FS−CNS2.

The ADCP measurements provided data about water depth and velocity. The observed water depth ranged from 21 to 68 m, with an increase of about 6–7 m from low to relatively high flow conditions. Table 4 lists the minimum, average, median and maximum value of the depth-averaged velocity along the ADCP transects as well as its standard deviation. No large variations were observed from low flow to relatively high flow conditions on minimum and average velocity, but the maximum depth-averaged velocity increased from 2.2 to 2.8 m/s. Figures 5 and 6 show the longitudinal distribution of the depth-averaged velocity in the ADCP transects collected in FS−CNS1 and FS−CNS2, respectively.

**Table 4.** Depth-averaged velocity in the ADCP transects collected in FS−CNS1 and FS−CNS2. Legend: *Vmin* = minimum depth-averaged velocity; *Vmean* = mean of the depth-averaged velocity; *Vst.dev.* = standard deviation of the depth-averaged velocity; *Vdepth-avg* = median of the depth-averaged velocity; *Vmax* = maximum depth-averaged velocity.


**Figure 5.** Depth-averaged velocity (blue) and bedline (black) for Long1\_11\_14\_000 transect (FS−CNS1).

**Figure 6.** Depth-averaged velocity (blue) and bedline (black) for Long1\_03\_05\_15\_000 transect (FS−CNS2).

Downstream from the mouth of the Solimões, the longitudinal transect was initially located on the rocky scour hole, with a depth larger than 60 m. After the scour hole, it was observed a wide sandstone terrace at depth of about 35 m, gently sloping downstream, followed by sharp 30 m decrease in depth. Downstream of this large depression, the bedforms were found [43]. The starting point of the bedforms is located within the flow recovery region, where the Amazon channel is widening, at approximately 4.7 km downstream of the confluence junction. Interestingly, the bedform size seems to increase and the shape changes with downstream distance.

Some parameters related to sediment transport were calculated. In alluvial channels, friction is related both to the grain resistance and to the form of bedform and the total shear stress is [2,45]:

$$
\pi\_b = \pi\_b' + \pi\_b'' \tag{1}
$$

where τ- *<sup>b</sup>* and τ*"b* are the skin friction shear stress and the form-related shear stress, respectively. *Water* **2020**, *12*, 1630

The skin friction shear stress was calculated as [46]:

$$
\pi\_b' = \rho \mathbb{C}\_d V\_{dep tl-av\chi}^2 \tag{2}
$$

where *Cd* is the drag coefficient, which was obtained as:

$$\mathbf{C}\_d = \frac{\kappa^2}{\ln^2(h/\varepsilon \,\mathbf{z}\_0)}\tag{3}$$

where κ is the von Kármán constant, *e* is the Euler number and *z*<sup>0</sup> is the zero−velocity height above the bed, which can be obtained as [47]:

$$z\_0 = 0.1 \, d\_{84} \tag{4}$$

where *d*<sup>84</sup> is bed grain diameter such that 84% of diameters are finer (Figure 7) [48]. The water density was calculated from the temperature *T* as:

$$
\rho = -0.0054 \ T^2 + 0.021 \ T + 1000 \ \text{(kg/m}^3\text{)}\tag{5}
$$

**Figure 7.** Grain size distribution [48].

The form-related shear stress was computed as [2]:

$$
\pi''\_{b} = \frac{1}{2} \rho V\_{depth-avg}^2 \frac{H\_{bf}^2}{L\_{bf}h} \tag{6}
$$

where *Hbf* and *Lbf* are bedforms height and length, respectively, and *h* is water depth.

Table 5 lists the minimum, average, standard deviation, median and maximum value of the bed shear stress along the ADCP transects calculated using Equation (1). No large variations were observed in the median value from low to relatively high flow conditions. Finally, the shear stress allows us to calculate the maximum suspended grain size *dss* [49]:

$$d\_{\rm ss} = \sqrt{\frac{18\rho\nu \cdot 0.8 \cdot \sqrt{\tau\_b/\rho}}{g(\rho\_s - \rho)}} \ (m) \tag{7}$$

where ρ*<sup>s</sup>* is the particle density, and ν is the water kinematic viscosity, that was calculated as:

$$\nu = \text{(5.85 } 10 - 10 \ast T^2) - \text{(4.85 } 10 - 8 \ast T + 1.74 \text{ 10} - 6) \text{ (m}^2\text{/s)}\tag{8}$$


**Table 5.** Bed shear stress τ*b*. Legend: τ*b-min* = minimum bed shear stress; τ*b-mean* = mean of the bed shear stress; τ*b-median* = median of the bed shear stress; τ*b-max* = maximum bed shear stress.

Table 6 lists the minimum, average, standard deviation, median and maximum value of the maximum suspended grain size along the ADCP transects calculated using Equation (7). The maximum suspended grain sizes were generally in the order of fine sand (0.125−0.250 mm) with a highest value in the range of medium sand (0.250−0.500 mm) in both flow conditions.

**Table 6.** Maximum suspended grain size *dss*. Legend: *dss-min* = minimum value of the maximum suspended grain size; *dss-mean* = mean of the maximum suspended grain size; *dss-st.dev* = standard deviation of the maximum suspended grain size; *dss-median* = median of the maximum suspended grain size; *dss-max* = maximum value of the maximum suspended grain size.


#### *3.2. Bedform Morphology. Results and Discussion*

Bedform characteristics were derived through several steps. The raw ADCP data were first extracted with WinRiver II. Then, the ADCP data were processed to get depth-averaged vertical and streamwise velocities according to the procedure described in Bahmanpouri et al. [50] with the addition of a further low-pass filtering to remove spikes and noise. The next stage was to track the bottom profile to detect the bedforms as anomaly relative to a reference depth. To this aim, a Matlab code was implemented: a seventh-grade polynomial fit curve was chosen (Figure 8) to define a reference depth used to detect using visual analysis individual bedforms as a succession of trough–crest–trough and to estimate their wavelength and wave height. Using this procedure, further metrics such as wave steepness, lee side and stoss side lengths and angles were calculated.

**Figure 8.** Polynomial seventh grade fit and longitudinal transects for (**a**) FS−CNS1 and (**b**) FS−CNS2.

Tables 7–9 list the minimum, average, median and maximum value of the wavelength λ*bf*, wave height *Hbf* and wave steepness *Hbf*/λ*bf* as well as their standard deviation for the bedforms observed during this study. The datasets for FS−CNS1 and FS−CNS2 were termed "Encontro das Aguas—1 November 2014" and "Encontro das Aguas—3 May 2015", respectively. Figures 9–11 show the frequency distribution for these parameters. For each parameter 14 classes in size were considered. The 2D analysis identified 36 and 70 bedforms, which were all classified from Table 1 [22] at least as *large dunes* (λ*bf* > 10 m), while *very large dunes* (λ*bf* > 100 m) were 86% and 70%, in FS−CNS1 and FS−CNS2, respectively. The minimum wavelength λ*bf* was equal to 55 m, observed in the FS−CNS2. The maximum wavelength was found in relatively high flow conditions and it was longer than 330 m. On average, the wavelength λ*bf* was 150 and 128 m for FS−CNS1 and FS−CNS2, respectively, while the average wave height *Hbf* was 3.7 m in both cases. The average steepness *Hbf*/λ*bf* was of 2.5% and 3.0% for FS−CNS1 and FS−CNS2, respectively. The distribution of the wavelength of the observed bedforms had a small percentage (14% and 30%) of bedforms shorter than 100 m.

**Table 7.** Wavelength λ*bf*. Legend: λ*bf*-*min* = minimum wavelength; λ*bf-mean* = mean of the wavelength; λ*bf-st.dev* = standard deviation of the wavelength; λ*bf-median* = median of the wavelength; λ*bf-max* = maximum wave length.


**Table 8.** Dune height *Hbf*. Legend: *Hbf-min* = minimum wave height; *Hbf-mean* = mean of the wave height; *Hbf-st.dev.* = standard deviation of the wave height; *Hbf-median* = median of the wave height; *Hbf-max* = maximum wave height.


**Table 9.** Dune steepness *Hbf*/λ*bf* (multiplied by 100). Legend: (*Hbf*/λ*bf*)*min* = minimum wave steepness; (*Hbf*/λ*bf*)*mean* = mean of the wave steepness; (*Hbf*/λ*bf*)*st.dev.* = standard deviation of the wave steepness; (*Hbf*/λ*bf*)*median* = median of the wave steepness; (*Hbf*/λ*bf*)*max* = maximum wave steepness.


**Figure 9.** Frequency distribution of bedform wavelength λ*bf* for (**a**) FS−CNS1 and (**b**) FS−CNS2.

**Figure 10.** Frequency distribution of bedform wave height *Hbf* for (**a**) FS−CNS1 and (**b**) FS−CNS2.

**Figure 11.** Frequency distribution of bedform steepness *Hbf*/λ*bf* for (**a**) FS−CNS1 and (**b**) FS−CNS2.

The lee side angle was constantly below 10◦, with a maximum value of 8.47◦ and 8.87◦ and an average value of 3.02◦ and 3.25◦, in FS−CNS1 and FS−CNS2, respectively. These values are generally consistent with those from large rivers where low-angle (<10◦) lee-side slopes are predominant [15]. The asymmetry, defined as the ratio of stoss side length to the bedform length, was on average 0.56 and 0.47 in FS−CNS1 and FS−CNS2, respectively.

In the rising stage from October 2014 (FS−CNS1) to April/May 2015 (FS−CNS2), on average, the wavelength decreased, the steepness increased and the wave height remained unchanged, while the maximum sizes increased. Furthermore, a comparison between the frequency distribution of bedform size in low and relatively high flow conditions showed an increase in wavelength and wave height as the river discharge increased, in agreement with the past literature studies (Figures 9–11). However, as the two ADCP transects are different in length of about 4 km and have a different number of bedforms, the comparison was repeated considering only the bedforms located on the same reach of the two longitudinal transects. The results confirmed the above findings. At the end, the bedforms

observed in this study were generally characterised by large wavelengths, ranging from 55 to 335 m, with a median value of 133 and 120 m, respectively, and wave height on average larger than 3 m. The wave steepness was in the range from 0.3% to 11%.

#### *3.3. Modeling Bedforms Morphology. A Comparison with Predictive Equations. Discussion*

As already pointed out, dune development is related to both their deformation during migration and their adaptation to flow variations [9]. Dune adaptation has been extensively investigated as an important process in river morphodynamics, but there is not yet a universal model to predict changes in dune sizes in response to flow variations. This morphological response has been often related to sediment mobility, which itself is a product of flow depth and velocity [9], and dune size has been related to flow depth as a result of interaction between large eddies in the flow and the sediment bed [17,51,52], Bedform wavelength was plotted versus bedform wave height [13,20,26,53] and the data were compared with the empirical relationships proposed by Flemming, which were based on 1491 deep sea, tidal and river bedforms [22]: 

$$H\_{bf} = 0.068 \,\lambda\_{bf}^{0.81} \,\mathrm{(m)}\tag{9}$$

$$H\_{bf-max} = 0.16 \,\lambda\_{bf}^{0.84} \,\mathrm{(m)}\tag{10}$$

where Equation (9) represents a range of steepness *Hbf*/λ*bf* from 0.08 to 0.1 [20]. Figure 12 shows the results for the bedforms observed in the field surveys, including their respective averages. The data were also compared with the equation proposed by Chen et al. [13] using experimental data collected in the middle–lower Changjiang (Yangtze) River (China):

$$H\_{bf} = 0.23 \,\lambda\_{bf}^{0.56} \, (\text{m}) \tag{11}$$

**Figure 12.** Bedform wavelength vs. bedform wave height.

Furthermore, the equation proposed by Lefebvre et al. [53], who used the data from the Rio Paranà (Argentina) [12] and those from the Lower Rhine (the Netherlands) [54], was included:

$$H\_{bf} = 0.13 \,\, \lambda\_{bf}^{0.59} \,\, (\text{m}) \tag{12}$$

Using the data from the two field surveys, it was possible to derive a new regression equation: ͷ

$$H\_{bf} = 0.21 \,\text{A}^{0.56}\_{bf} \,\text{[m]} \tag{13}$$

which is very close to Equation (11). Furthermore, theory and laboratory studies in uniform and steady flow suggest that dunes with a steepness *Hbf*/λ*bf* less than 0.06 are either non-equilibrium bedforms or represent an equilibrium adjustment of the bed form, in which maximum steepness is precluded by hydraulic constraints, notably a depth limitation [20]. The equilibrium line corresponding to *Hbf*/λ*bf* = 0.06 was also included in Figure 12. All the bedforms, but one, were below Flemming's maximum line (Equation (10)). Only some of *large dunes* (10.0 > λ*bf* > 100 m) were well aligned between Equation (11) and Equation (9) and close to the equilibrium line, while in many cases they showed a large scatter from these lines. On the other hand, most of the *very large dunes* (λ*bf* > 100 m) were well aligned with both Chen et al. [13] (Equation (11)) and Flemming's (Equation (9)) lines and close to the equilibrium line, but several bedforms from FS−CNS2 and also a number from FS−CNS1 had a low wave height, corresponding to a steepness in the order of 0.01–0.02, so they may represent bedforms in adaptation. It is worth noting that during FS-CNS1, the longitudinal transect was collected after seven days of near-constant low-flow discharges, while FS-CNS2 was conducted during a period of continuously rising flow discharges. This could explain why the dune field may have mostly obtained a stable equilibrium with the flow conditions during FS-CNS1, while dune field was in a transitional phase during FS-CNS2 as it adjusted to the increasing flow discharge.

The length of the bedforms observed in the three field surveys were plotted against their steepness *Hbf*/λ*bf* in Figure 13 and compared with the relationship proposed by Carling et al. [20]:

$$\frac{H\_{bf}}{\lambda\_{bf}} = 0.1027 \,\lambda\_{bf}^{-0.615} \, (-) \tag{14}$$

and a new regression equation was derived:

$$\frac{H\_{bf}}{\lambda\_{bf}} = 0.208 \ \lambda\_{bf}^{-0.437} \ (-) \tag{15}$$

**Figure 13.** Bedform wavelength vs. bedform steepness.

Several bedforms, in both flow conditions, were characterised from steepness in the order of 0.01–0.2. This indicates that the bedforms with low steepness observed in FS−CNS2 were developing with the increasing discharge, while, on the other side, those in FS−CNS1 were in the process of crest flattening and elongation, having been formed during the previous high flow conditions [13].

As already mentioned above, dune sizes are often thought to scale with flow depth [10,11,17,26,51,52,55], as developing dunes cannot emerge out of the water [26], but, in the literature, other scaling relationships based on depth and grain size [24], transport stage [23,56] and transport stage and Froude number [57] were proposed. Transport stage is generally defined as any metric that is composed of a ratio of the shear stress to a grain size [25,26], including the Shields number and the Rouse number, which is defined as:

$$Ro = \frac{w\_{\\$}}{\text{\textsuperscript{K}\,\,\,u^{\*}}} \tag{16}$$

where *ws* is particle settling velocity and *u\** is the shear velocity, which can be obtained from the total bed shear stress τ*<sup>b</sup>* as *u\** = (τ*b*/ρ) 0.5. The most widely applied scaling equation for dune size is that of Yalin [51], where bedform wavelength and wave height are related to the flow depth as:

$$
\lambda\_{bf} = 5 \,\text{h (m)}\tag{17}
$$

$$H\_{bf} = \frac{h}{6} \text{ (m)}\tag{18}$$

while Yalin [52] suggested a theoretical value of λ*bf* = 5 *h* for equilibrium dunes in deep flows. Venditti [55] analyzed Allen's [58] dataset and identified a range of variability between *h* and 16 *h* for bedform length and between 1/40 *h* and 1/6 *h* for the bedform height. He also argued that the reason for this variability is that bedform sizes are dependent on transport stage and lag changes in the flow conditions.

Bradley and Venditti [26] re-evaluated seven predictive equations, including those from Yalin, linking dune dimensions to other variables such as flow depth, grain size, transport stage and Froude number. The data compilation using 498 observations coming from 21 flume experiments and 20 field studies shows that dune height and length follow a power law:

$$H\_{bf} = 0.051 \,\, \lambda\_{bf}^{0.77} \,\, (\text{m}) \tag{19}$$

which is very similar to Flemming's equation (Equation (9)). Most of the data for bedform length and height were ranging from *h* to 16 *h* and from *h*/20 to *h*/2.5, respectively. Bradley and Venditti [26] found that the predictive power of all the scaling relations was generally poor probably because these relations are not able to capture any effect of the flow variability over the time. They also observed that dunes in smaller channels conform to a different height scaling than dunes in larger channels, which reflects a change in dune morphology from strongly asymmetric dunes with high lee angles in flows <2.5 m deep to more symmetric, lower lee angle dunes in flows >2.5 m deep [26]. This implies a different process control rather than a continuum of processes as depth increases [26]. Hence, from the analysis of only data in channels with depth *h* > 2.5 m, they derived two different equations for wave height, a linear regression equation and a non-parametric scaling equation:

$$H\_{bf} = 0.13 \, h^{0.94} \text{ (m)}\tag{20}$$

$$H\_{bf} = \frac{h}{7.7} \text{ (m)}\tag{21}$$

while, for all the data, they derived for wavelength again a linear regression equation and a non-parametric scaling equation:

$$
\lambda\_{bf} = 5.22 \, h^{0.95} \, (\text{m}) \tag{22}
$$

$$
\lambda\_{bf} = 5.9 \, h \, (\text{m}) \tag{23}
$$

Bradley and Venditti [26] recommended to apply non-parametric scaling relations, i.e., Equations (21) and (23), concluding that if it is clear that dunes increase its sizes with the scale of the river system, a sound explanation for how flow depth rules the equilibrium sizes of dunes is still lacking. Rather, the apparent scaling of dunes with flow depth may be only indirect and emerge because shear stress/velocity, both depending on depth, play a key role in dune morphology [26].

Figures 14 and 15 present the relationship between water depth and bedform wavelength and wave height, respectively, for the dunes observed in this study. The data were not exceeding the upper limits of scaling reported by [55]. On the other hand, while length data were quite well aligned between the curves for λ*bf* = 5.9 *h* (Equation (23)) and λ*bf* = *h* (Figure 14), height data ranged mostly from *Hbf* = *h*/6 to *Hbf* = *h*/40, with some values above and below these curves (Figure 15). The bedforms with low wave height, i.e., *Hbf* < *h*/20 or *Hbf* < *h*/40, observed mostly in FS−CNS2 were probably adapting to the raising stage to high flow conditions. Finally, most of the wave height data were below the curve for *h*/10, confirming the findings observed in large rivers where height is often only 10% of the local flow depth [15].

**Figure 14.** Bedform wavelength vs. water depth.

**Figure 15.** Bedform wave heigth vs. water depth.

Bedform sizes were also related to the transport stage expressed by the suspension number *u\**/*ws* [25], which is the inverse of the Rouse number (Equation (15)) multiplied by the von Kármán constant. Figure 16 presents the distribution of the relative bedform wavelength λ*bf*/*h* and wave height *Hbf*/*h* against the suspension number *u\**/*ws*, being in both flow conditions, in the range from 0.3 to 0.6. The lee side angle was below 10◦, with an average value of 3.02◦ and 3.25◦ in 2014 and 2015 surveys, respectively. Lee side angles showed an interesting relationship to wave steepness (Figure 17). Bedform steepness grew gently with lee side angle and became constant above 6◦ at *Hbf*/λ*bf* = 0.05, suggesting an interrelation between these parameters in this case [53]. Finally, the analysis of vertical velocity measured in the ADCP transects showed no flow separation due to the low lee side angles, as reported in past studies on low-angle dunes [1,15,59].

**Figure 16.** Suspension number *u\**/*ws* vs. λ*bf*/*h* (left) and vs. *Hbf*/*h* (right).

**Figure 17.** Bedform lee face angle vs. bedform steepness.

#### *3.4. Comparison with Literature Data Sets from Large Rivers. Discussion*

The field data collected is this study were compared with those from some literature data sets collected in large rivers. These data sets are: the side-scan sonar and sub-bottom profiler data collected in August 2003 in the middle–lower Changjiang (Yangtze) River (China) [13], and the multibeam echosounder data collected in January 2004 in the Lower Rhine (The Netherlands) [54], in May 2004 in the Rio Paraná (Argentina) [12], and in April 2008 in the Empire Reach of the lowermost Mississippi River (USA) [60]. Table 10 lists the main size parameters for these data sets.

**Table 10.** Main parameters of bedforms from the literature data sets. Legend: λ*bf-mean* = mean of the wavelength; *Hbf-mean* = mean of the wave height; (*Hbf*/λ*bf*)*mean* = mean of the wave steepness.


Figure 18 presents the relationship between bedform wavelength and bedform wave height for the dunes observed in the present study and those reported in the above literature data sets along with the empirical regression laws in Equations (9)–(13). In the range of *medium dunes* (5.0 > λ*bf* > 10 m) or more, the data from the Lower Rhine and Mississippi River were aligned with Equation (9), while the data from the Changjiang River were proximal to Flemming's maximum line (Equation (10)) and their steepness was larger than 0.06. Moreover, several data from the Lower Rhine were below Flemming's line (Equation (9)). In the range of *large dunes* (10.0 > λ*bf* > 100 m), most of the dunes from Rio Paraná and several from the Changjiang River were aligned with Equation (13) and close to the equilibrium line, but a number of dunes from Rio Paraná were below Flemming's line (Equation (9)). In the range of *very large dunes* (λ*bf* > 100 m), most of the literature data were close to the Amazon River data and it is possible to observe that a number of data from the Changjiang River were characterised from low wave height and steepness in the order of 0.01, as they were probably adapting to changes in flow conditions.

**Figure 18.** Bedform wavelength vs. bedform wave height, comparison with literature data from large rivers.

#### **4. Conclusions**

In alluvial rivers, the geometry of the bed topography is the result of a complex interaction among several hydrodynamics and sedimentary processes acting under the constraint of varying boundary conditions. In sand-bedded alluvial channels, the bottom boundary consists of a labile bed comprising bedforms of many different scales and geometries. Bedforms are deformations of a sand bed that are smaller than channel-scale bar forms and that have specific geometric properties [26,55]. They also are subjected to deformation and adaptation to changes in river flow [9]. This paper presented the results of a study about the morphology of the bedforms observed in the area of the Negro/Solimões confluence. Two surveys using acoustic Doppler velocity profiling (ADCP) were carried out in low flow (FS-CNS1) conditions, after seven days of near-constant discharge, and in relatively high flow (FS-CN2) conditions, during a period of continuously rising flow discharges. The observed bedforms were mostly in the range of *large* and *very large dunes* according to Ashley [22] classification with a maximum wavelength and wave height of 350 and 12 m, respectively. Second, during (FS-CN2), maximum bedform sizes as well as in the frequency distribution of bedform size were comparatively larger, as could be expected from the literature studies. Third, most of the *large dunes* (10.0 > λ*bf* > 100 m) and *very large dunes* (λ*bf* > 100 m) were generally in equilibrium with flow conditions. On the other side, some bedforms observed in relatively high flow conditions were developing to adjust to the continuously increasing flow discharge, while only some in low flow conditions were in the process of crest flattening and elongation, having been formed during the previous high flow conditions. Fourth, the data were within the upper limits of scaling with water depth reported in the literature, but, while length data were quite well aligned with scaling curves, height data showed a scatter from these curves. Fifth, in both surveys, the lee side angle was below 10◦ with an average value of about 3.0◦ and no flow separation was observed confirming recent literature studies on large rivers [15], and wave steepness grew gently with lee side angle and became constant above 6◦ at *Hbf*/λ*bf* = 0.05. Finally, a comparison between the data collected in this study and past literature studies on large rivers suggested that a number of the largest bedforms were probably adapting to discharge variations in the river.

**Author Contributions:** All the authors contributed to this study. Analysis and investigation, C.G., I.M., N.P.F.J. and M.I.; data processing, C.G. and I.M.; writing—original draft preparation, C.G.; writing—review and editing, I.M., N.P.F.J. and M.I.; visualization, C.G., I.M. and M.I.; funding acquisition, C.G., N.P.F.J. and M.I. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was carried out within the Clim-Amazon Research Project funded by grant agreement FP7 INCO-LAB n◦295091 from the European Commission. Naziano Filizola acknowledges funding by CAPES-Procad Amazonia and CNPq-AmazonGeoSed Project.

**Acknowledgments:** The authors acknowledge Mark Trevethan for his valuable collaboration in the Clim-Amazon Project and the CPRM (Geological Survey of Brasil) for supplying the research vessel, instrumentation, technical assistance with sampling; Bosco Alfenas, Andrè Martinelli Santos, Daniel Moreira, Arthur Pinheiro, Paulo Melo, Nilda Pantoja, Andre Zumak and João Andrade for their assistance with sampling during the surveys; and finally Farhad Bahmanpouri for his assistance in the extraction of the ADCP data.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
