Scanning PIV of turbulent flows over and through rough porous beds using refractive index matching

This paper presents image velocimetry measurements on turbulent flows adjacent to a permeable bed made of randomly packed glass particles. For measuring flow velocities inside the bed, the refractive index of the glass particles was matched with that of the fluid. By continuously scanning in the transverse direction, we measured the streamwise and vertical velocity components within a three-dimensional domain (3D2C-PIV), including firstand second-order turbulent statistics. We established how the scanning travel speed is associated with the laser sheet thickness and the space-time velocity fluctuations for collecting reliable measurements. The methodology was applied to free-surface flows over a sloping bed under low relative submergence and supercritical conditions. Spaceand time-averaged profiles were obtained in a representative elementary volume as defined by the double-averaging procedure (Nikora et al. in J Hydraulic Eng.127(2):123–133, 2001). A turbulent boundary layer over the rough bed was observed when experiments were run at intermediate Reynolds numbers Re = O(1000) . Apart from measuring subsurface velocities, this method shed light on the part played by the rough bed in the overall flow dynamics: the roughness layer was a buffer region within which porosity varied sharply and turbulent stress was rapidly dampened. * Gauthier Rousseau gauthier.rousseau@gmail.com 1 Swiss Federal Institute of Technology in Lausanne, EPFL, Laboratory of Environmental Hydraulics, EPFL-ENACLHE, Station 18, 1015 Lausanne, Switzerland Experiments in Fluids (2020) 61:172 1 3 172 Page 2 of 24


Introduction
Knowledge about how boundary-layer flows are affected by porous boundaries is central to many fluid applications, ranging from shear flows of air over and through forest canopies or foams, to heat transfer problems (Ghisalberti and Nepf 2006;Mahjoob and Vafai 2008;Suga et al. 2010). Open-channel flows over rough permeable beds are another case in point. Although these flows have been extensively studied (Nezu and Rodi 1986;Nezu 2005), there remain some unanswered questions. For instance, mountain rivers exhibit two distinctive features not shared by lowland rivers: (i) their flow-depth and roughness scales are similar, which makes their turbulent structures far more complex, and (ii) their bed permeability is much higher. Both features help to explain why classic boundary layer assumptions exhibit strong limitations for predicting river processes such as flow resistance, hyporheic exchanges or the incipient motion of sediment (Recking 2009;Rickenmann and Recking 2011;Boano et al. 2014;Prancevic and Lamb, 2015;Lamb et al. 2017;Rousseau 2019).
Making a grain-scale examination of the physical processes involved (see Fig. 1) reveals how bed protuberances act as obstacles to flow: they create local wakes, promote vorticity and produce different types of turbulent structures. As a result, turbulent boundary layers show large spatial heterogeneity (Mignot et al. 2009a). Flow paths-whether in the roughness or subsurface layer-are tortuous. (These paths are sketched as dashed arrows in Fig. 1.) All these elements make it very difficult to use classic eddy-resolving modelling methods at the river scale essentially owing to their high computational cost (Keylock 2015).
A fine-grained description of turbulent flows adjacent to complex porous boundaries seems currently out of reach. However, this does not mean that we cannot gain any theoretical insight into the issue. Working at the mesoscopic scale, averaging flow properties over a representative elementary volume and taking their time averages make it possible to provide a consistent physical picture. This approach relies on the double-averaging concept developed by Nikora et al. (2001Nikora et al. ( , 2007 for studying open-channel flows over rough beds: an approach that has taken inspiration from previous contributions examining flow over vegetated canopies (Wilson and Shaw 1977). For these flows, bed characteristics vary along the vertical axis and flow fluctuates in time and space making the problem more complicated to treat than for homogeneous materials. By averaging the Navier-Stokes equations over time, then averaging the resulting equations over a thin volume parallel to the bed surface, one obtains the double-averaged momentum equations. This procedure is reminiscent of the Reynolds decomposition used for obtaining the Reynolds-averaged Navier-Stokes equations (RANS). As with the RANS equations, time-averaging produces the Reynolds tensor which is interpreted as the turbulent stress tensor. Additional terms arise when space-averaging the momentum balance equations: the viscous drag, the pressure drag and the dispersive stress (also termed the form induced stress). Dispersive and turbulent stresses are algebraically defined, and thus, they can be determined from experiments or eddy resolving numerical simulations. Dispersive stresses are associated with the spatial variability of the velocity field, and their study is more delicate. Recent investigations have revealed that their contribution to the momentum balance might have a substantial effect at the bed interface (Voermans et al. 2017;Fang et al. 2018).
As shown in Fig. 1, spatial averaging defines a porosity profile, which can be used to distinguish between three specific regions: the surface layer, the roughness layer and the subsurface layer. Although porosity profiles are continuous for natural rough beds, the technical literature devoted to flows adjacent to permeable walls generally considers a discontinuity, i.e. a porosity jump from a finite value to unity at the bed interface (Beavers and Joseph 1967;Mendoza and Zhou 1992;Breugem et al. 2006;Tilton and Cortelezzi 2008;Rosti et al. 2015;Zampogna and Bottaro 2016;Lācis and Bagheri 2017). This porosity jump is associated with a momentum transfer, whose strength can be estimated by considering a Brinkman correction, that is, a model describing how viscous stress propagates across the layer with a jump condition at the interface (Brinkman 1949;Ochoa-Tapia and Whitaker 1995). When the flow is turbulent in the close vicinity of the interface, the flow structure close to the bed is substantially altered by permeability. For instance, turbulent eddy vortices may penetrate into the bed (Suga et al. 2010). This scenario is ubiquitous in nature (flows in rivers, over canopies or adjacent to biological surfaces such as feathers, for instance). Recent contributions demonstrated that all these flows exhibit a similar behaviour (Ghisalberti 2009;Keylock et al. 2019;Bottaro 2019;Chagot et al. 2020). However, there is no consensus for modelling such flows and this starts by notifying the existence of at least two school of thoughts for their theoretical examination: the homogenisation theory (Bottaro 2019) and the double averaging methodology (Nikora et al. 2007). This lack of consensus is likely to reflect the dearth of high-resolution experimental data ( e.g. see discussion in Cameron et al. (2017)- § 2.1 ), which results from the difficulties of probing velocities over and through permeable beds without disturbing flows. Noninvasive methods such as Particle Image Velocimetry (PIV) or Laser Doppler Velocimetry (LDV) are able to bring quality experimental data over the permeable bed (Mignot et al. 2009a;Suga et al. 2010;Manes et al. 2011;Blois et al. 2014;Cameron et al. 2017;Efstathiou and Luhar 2018), but are not adapted to capture flows through opaque beds. Authors such as Pokrajac and Manes (2009) ;Suga et al. (2018); and Grain-scale processes of a turbulent flow over a rough permeable bed and double-averaged quantities. The porosity and velocity profiles result from the double averaging procedure over a thin layer parallel to the mean bed surface over the length L. The flow is subdivided into three specific regions: the surface layer, the roughness layer and the subsurface layer. The roughness layer is bounded by: (i) the roughness crest z rc above which the averaged porosity is unity, and (ii) the troughs of the roughness elements z t , where the bed porosity b is reached. The red dotted arrows represent streamlines through the roughness and subsurface layers forming the permeable bed Rouzes et al. (2019) get around this issue by creating artificial beds that enabled the visualisation of flows through the flume's sidewall, but these bed structures are far removed from real-world beds.
In order to capture flows adjacent to permeable beds exhibiting similarities with real-world scenarios, this paper presents an experimental procedure for investigating open channel flows over and through a randomly packed bed of spherical glass beads under low relative submergence conditions, i.e. with a bed roughness that is comparable in size to flow depth. Further physical insights into flow behaviour are provided by using the double averaging concept. This approach focuses on determining the porosity and velocity profiles as well as the dispersive and turbulent stresses present at the mesoscopic scale in a representative elementary volume, that is, whose typical length is about ten particle diameters. Our protocol is innovative in that it collects flow and material characteristics by coupling the Refractive Index-Matched Scanning (RIMS) method with a PIV technique to measure flow velocities. Refractive index matching (RIM) involves matching the fluid's and the beads' refractive indices so that the mixture becomes transparent (see Fig. 2). Several authors have used RIM for visualising flows: for instance, Hassan and Dominguez-Ontiveros (2008) used it for measuring interstitial velocities in pebble reactors, while Aussillous et al. (2013) determined grain velocities in laminar flows. Budwig (1994) and Wiederseiner et al. (2011) reviewed the different RIM techniques. Refractive Index-Matched Scanning (RIMS) combines RIM and a scanning procedure: a laser sheet is displaced through the index-matched medium to collect information from successive measurement planes. Most RIMS-based methods have been used as tomography, i.e. a technique for locating solid elements in a three-dimensional domain (Huang et al. 2008;Dijksman et al. 2012;van der Vaart et al. 2015).
With similar objectives to ours, Voermans et al. (2017) used an index matching technique to study mass and momentum transfers across the roughness layer (termed sediment-water interface). Voermans et al. (2017) produced porosity and double averaged profiles from several fixed measurement planes. As they did not provide a detailed experimental protocol, there is little information on the measurement accuracy, in particular, reproducibility. Ni and Capart (2015) simultaneously scanned and measured velocities in a saturated granular flow (erosion of a granular bed by a turbulent flow). Measuring solid and liquid displacements was difficult in such a flow, because they had to be captured at the same time. With the experimental set-up presented here, we do not face the same difficulty because the granular bed remained static.
Reviewing the technical literature with a focus on methods coupling PIV and scanning shows the predominance of sophisticated procedures involving rotating scanners (Brücker 1997) or stereoscopic PIV with a monitored rotating mirror (Hori and Sakakibara 2004). In these scanning experiments, authors determined velocities from consecutive images taken in the same measurement plane, and they tried to deduce the instantaneous three-dimensional (3D) velocity field by quickly moving the laser sheet, which required high-speed cameras (operated at more than 2000 frames per second) and specific equipments. The approach presented in this paper was different, less complex and less costly in comparison. Our goal was to determine time-and space-averaged quantities that were compatible with a certain range of scanning speed. Scans were performed at constant speed, which means that two consecutive measurement planes were not exactly at the same transverse coordinate. This procedure turned out to be convenient, but admittedly, it is unconventional within the PIV community. As a consequence, part of the paper is devoted to explaining how double-averaged profiles can be collected by scanning the medium continuously and thereby provide a wealth of information on the flow turbulence characteristics. The procedure enabled us to measure the streamwise and vertical velocity component in Fig. 2 Left: Two beakers containing equal quantities of borosilicate beads representing the porous bed. The left-hand beaker contains water and the refractive index mismatch means that the beads are visible; the right-hand beaker contains a fluid whose refractive index matches that of the beads, rendering them invisible, but the interstitial fluid can be probed using flow visualisation techniques. Right: Pho-tograph of a gravity-driven flow (on a i = 0.5% slope) over a porous bed made of the same borosilicate beads. The RIMS technique enables us to visualise the interior of the roughness and subsurface layers: the black disks are the borosilicate beads illuminated by the laser sheet, whereas the small dots are tracers. A PIV technique was used to measure the flow velocity a three-dimensional domain (a method often termed 3D2C-PIV). To the best of our knowledge, coupling simultaneously PIV and transverse scanning has not so far been presented for studying flow/porous structures interactions. While the procedure does not necessarily reduce the amount of data to be handled, we see it as a fast and convenient method, especially when experimental conditions impose a short time interval for performing measurements.
The article is organised as follows: Sect. 2 describes the experimental set-up and protocol used in our experiments. Section 3 explains the scanning procedure used to determine flow properties at the mesoscopic scale. We show that this procedure imposed constraints on the transverse scanning speed relative to the laser sheet thickness, spatial variability and turbulence. An empirical validation is made by comparing the results from a fixed measurement plane and the data obtained from the same plane, but with the PIV-RIMS procedure. Section 4 focuses on flow uniformity, data reproducibility and uncertainties. Section 5 presents preliminary observations based on the PIV-RIMS method which allows capturing flows from the subsurface to the free surface.

Flume and materials
The experiments were performed in a 6-cm wide, 2.5-m long flume with an adjustable slope i, as shown in Fig. 3. A constant head tank provided a steady fluid discharge into the system. Equal proportions in mass of borosilicate beads of two diameters (7 and 9 mm) were randomly packed into the flume bottom, forming the coarse-grained bed. Median particle diameter was thus d p = 8 mm . A bed composed of beads of the same diameter would arrange itself in parallel layers causing undesirable bias in the averaged porosity and velocity profiles. Before each run, the upper layer was flattened out to form a uniform bed of height h s = 5 cm . Flow disturbances at the flume inlet were reduced using straighteners, and the region of interest (ROI-where measurements were taken) was located far upstream of the permeable grid placed at the flume outlet to maintain the beads while letting the flow seep out of the bed. Even though surface flows were supercritical, the downstream condition at the flume outlet affected the flow dynamics. For example, if the grid had been replaced with an impermeable wall, a dead zone would have appeared in the bed just upstream of the wall, causing flow resurgence and changes to the surface flow. Despite this measure, we could not exclude the development of substantial pore-pressure variations near the flume outlet. This problem is discussed in Sect. 4.3.
The iso-index fluid was prepared by mixing volumetric concentrations of 40% ethanol and 60% benzyl alcohol. The refractive index n f of the resulting fluid matched that of the borosilicate beads. Using a digital refractometer (ATAGO RX-5000 ), we found n f = n borosilicate−glass = 1.472 ± 0.002 at 20 • C . The iso-index fluid's physicochemical characteristics were close to those of water. Using a Cannon-Ubbelhode viscometer, we measured its kinematic viscosity at a temperature of 20 • C : f = 3.0 ± 0.1 mP ⋅ s . Its density was f = 950 ± 10 kg ⋅ m −3 . (Details of these measurements are provided in Rousseau (2019)-Annex D.) These values were close to those obtained by Chen et al. (2012). According to these authors, surface tension was about f = 31 ± 1 mN ⋅ m −1 , which is a factor of 2 lower than that of water. Surface tension was thus assumed to have negligible effects on our experimental flow dynamics. Fluid and sediment characteristics are measured at 20 • C and are summarised in Table 1.
The borosilicate beads' density was s = 2200 kg m −3 . Compared to other materials used in RIM techniques, combining borosilicate, ethanol and benzyl alcohol leads to mixtures whose relative density is close to that found in real-world scenarios like river engineering. Meeting these similarity criteria (e.g. the Shields and Froude numbers) is necessary to obtain flow conditions that mimic those encountered in real-world scenarios (e.g. a shallow flow on a steep slope with a low sediment transport rate, such as mountain streams). If we had followed (Ni and Capart 2015) and used a sediment of poly(methyl methacrylate), with a density of s = 1190 kg ⋅ m −3 , it would have been impossible to conduct experiments on steep slopes without observing sediment transport. The same observations would be expected when employing a popular RIM fluid mixture made of NaI (e.g. as in Voermans et al. (2017)) since, in this case, the fluid has an unexpectedly high density: f ,NaI = 1770 ± 10 kg ⋅ m −3 . In the present context, we needed a stable bed which could resist the stream's erosive action when the flow reached supercritical states.
The iso-index fluid was initially contained in a reservoir (with a volume of about 20 L) connected to a second reservoir below whose level was maintained constant using an overflow pipe, ensuring a steady flow rate into the flume with a flow depth in the centimetric range. Two valves controlled the desired flow rate: the first valve was manually controlled and regulated the base flow. The second was driven by an electro-valve and was used to adjust the flow rate to the desired value Q f . As shown in Fig. 3a, the reservoir was fixed at the flume's upstream end to obtain constant pressure heads regardless of the flume inclination. As the inclination did not exceed 8%, it had a negligible influence on the pressure head. Uncertainties on the flow rate were lower than 5%. The total discharge Q f was about 200 mL/s. With a limited volume of fluid per run, the steady state was reached during a small time interval of about 40 s. This short time interval pushed us to find an alternative to the usual fixed-laser-sheet technique by developing the PIV-RIMS protocol.
The iso-index fluid was chemically stable. At the interface between the flow and air, ethanol evaporated, a problem which might affect the fluid's refractive index in the long run. Ethanol was thus added between two consecutive runs. The fluid mixture was mixed with a stick, and the refractive index n f controlled to stay at the matching value. It must be emphasised for the reader interested in this mixture that the set-up requires specific security equipments (ventilation hood) because of ethanol evaporation and benzyl-alcool dissolving capacity. Small quantities of a fluorescent dye (Rhodamine B) were added to the fluid to increase the contrast between the beads and fluid. Our laboratory had previously used this combination of borosilicate beads and Rhodamine B to determine bead positions in tomography experiments on particle segregation in granular flows (van der Vaart et al. 2015).

Optical system
Frame sequences were recorded using a Basler acA2040-180kc camera operated at a rate of 420 frames per second and a resolution of 1496 × 700 pixels (px), giving an inter-framing time of Δ t = 2.3 ms. The lens' focal length was 35 mm, and the aperture was f/2.8. The camera was placed 30 cm from the sidewall, giving a field of vision of 73.8 × 34.5 mm 2 . Thus, the mesoscopic scale L over which spatial averaging was performed was about 8 cm or ten bead diameters. The flow was seeded with micrometric PIV tracers (hollow borosilicate glass spheres 8 − 12μm in diameter with a density of 1.1 10 3 kg∕m 3 ). We observed no particle deposition during the runs, a phenomenon sometimes reported in other RIM experiments (Voermans et al. 2017). This is probably because seeding particles and beads were both made of borosilicate glass. The tracers were lit up by a 4-W diode-pumped solidstate laser (emitting at 532 nm) mounted on a linear unit. The laser sheet's linear movement perpendicular to the flow enabled us to scan the ROI by taking images from the flume's sidewall. The laser sheet thickness was estimated at w LST = 1 − 2 mm. Measurements were taken in different transverse planes (along the y-axis) using a motorised carriage (see Fig. 3).
Here, we provide an estimate of the minimal length scale that can be resolved in our experimental set-up. We follow the procedure proposed by Kähler et al. (2016, pp. 18-19). An interrogation window of 32 × 32 px has a size of approximately 1.6 × 1.6 mm 2 in the physical space. By assuming typical scales of surface velocities ū ∼ 0.5 m s −1 and flow depth H = 1 cm, we deduce the rate of energy dissipation per mass unit: r =ū 3 ∕H ∼ 12.5 m 2 ∕s 3 . This leads to a Kolmogorov length scale of = ( 3 f ∕r) 1∕4 ∼ 40 m and a Taylor scale of g = √ 10 2∕3 L 1∕3 ∼ 0.7 mm. The resolution does not allow us to resolve the Kolmogorov scale, a situation quite common in PIV measurements (Kähler et al. 2016). Spatial resolution is on the order of the Taylor length, which implies that most vortices in the turbulent spectrum can be resolved.

Transverse scanning, porosity profiles and the vertical origin
Shifting the laser along the y-axis made it possible to take images in parallel planes and thus infer bead positions (x b , y b , z b ) n and diameters D n . After determining the bead positions, we built up a three-dimensional matrix of the porosity for the ROI, with a resolution of approximately one tenth of a bead diameter. This porosity array generalised the roughness geometry function defined by Nikora et al. (2001)). Each entry took the value of 0 if the datapoint (x, y, z) lay in a bead and 1 if it did not. When the point was close to the bead-fluid interface, the entry took a value ranging from 0 to 1 representing the volume-averaged porosity of the cell centred at (x, y, z) . We then obtained the averaged porosity for a slice of the porous bed at a position y on the laser sheet by summing the array over i and then dividing by the window length. The discrete spatial averaging was defined by: Similarly, we defined a cross-stream-averaged porosity profile by averaging B in the x-and y-directions (see Fig. 4): (1)

Fig. 4
Porosity measurement using RIMS: once the beads have been located, we define the porosity array B(x, y, z). a The B field averaged in the x-direction in a slice located at Y m = y − y w = 25 mm from the wall position y w . b Averaging B over x and y gives us a smooth porosity profile. The relative vertical coordinate is denoted by z � =0.8 = z − z b and is computed from z b , which is the vertical coordinate for which bed porosity is 0.8. The flow depth was defined by h f = z surf − z b , where z surf is the free-surface position and z b is the bed level. z surf is determined by extracting the surface position from each frame and then by performing an averaging over time. For low relative submergence flow conditions, i.e. for S m = h f ∕d p ∼ 1 , a definition of z b is essential. This quantity is used to deduce important flow parameters such as flow depth h f and, in consequence, the bed shear stress ( b = f gh f sin( ) ). A slight change in the definition may significantly alter these flow characteristics and the interpretation of the results, as highlighted by Pokrajac et al. (2006).
Here, z b is given by z =0.8 , which is the vertical coordinate where porosity is equal to 0.8. This gave a z b slightly below the roughness crest z rc . As there is no consensus about the definition of z b for rough beds, this choice might seem arbitrary and unconventional since z b is commonly given by the roughness crest z rc (e.g. Pokrajac et al. (2006); Cameron et al. (2017); Fang et al. (2018)). However, z =0.8 had the advantage of providing consistent comparisons between profiles at the mesoscopic scale when random packed bed arrangement was changed. We take a closer look at the issue of this choice in Sect. 4.2. In the following sections, the vertical coordinate is always referred to as the bed level, i.e.
Another possibility for the vertical origin found in the literature is the zero-plane displacement, which is commonly deduced by fitting the differential form of the log velocity profile to the experimental profile (e.g. Cameron et al. (2017)). We did not choose this option because a log profile is not expected to occur under low submergence conditions or when free surface level is below the roughness crest as commonly observed in gravel-bed rivers (Jiménez 2004).

Image velocimetry processing
The selection of an appropriate image processing technique for our experiments was constrained by two factors: (i) we estimated that the fluid passing through our roughness layer would exhibit substantial velocity variations, and (ii) we knew the gaps between the beads were narrow. These features would require the use of image velocimetry tools able to function across a sizeable dynamic range. We tested different methods, from classic PIV to the more elaborate particle tracking velocimetry (PTV). The open-source Python library, openCV, was the best suited to our needs.
The algorithm measures the local optical flow by means of a pyramidal application of the Lukas-Kanade method (Bouguet 2001). The optical flow method obtains the displacement field by minimising the squared Displaced Frame Difference (DFD). The methodology is similar to that of PIV algorithms, but it is optimised to be able to extract the displacement of any feature. Indeed, in classic PIV, algorithms are optimised for measuring particle displacement. For a better sense of the equations underpinning the algorithm, and its difference from classic PIV, the reader is referred to Liu and Shen (2008); Heitz et al. (2010);Boutier (2012). In turbulence, this methodology was used by Miozzi et al. (2008) and more recently by Zhang and Chanson (2018). Appendix 1 provides further information on the image velocimetry techniques used in our experiments. This appendix includes a test on the Case-A proposed in the 4th PIV challenge (Kähler et al. 2016) which presented similar challenges to ours, i.e. large velocity gradients and high dynamic velocity range. The algorithms showing the test results are all available from the public GitHub repository: https ://githu b.com/grous sea/opyfl ow.

Quantities of interest within the double-averaging framework
Turbulent flows over rough permeable beds exhibit strong spatial and temporal variability. The double-averaging concept was developed to cope with flow variability (Nikora et al. 2007). We consider a steady uniform flow and seek to define its mesoscopic flow properties.
Step 1 involves using the generalised Reynolds decomposition by breaking down the local instantaneous velocity into the time-space-averaged value ⟨u k ⟩ (with k = x , y or z), the local disturbance ũ k and the temporal fluctuations u ′ k in the three spatial directions. For a two-dimensional open-channel flow, the double-averaged decomposition gives: The superscript • represents time averaging, the brackets ⟨•⟩ are the intrinsic space averaging (i.e. over the fluid phase only), and the tilde superscript • denotes the local spatial disturbance. The control volume's dimensions in the x-and y-directions are sufficiently large for the mean fluctuating velocities ⟨ũ x ⟩ , ⟨ũ y ⟩ and ⟨ũ z ⟩ to be negligibly small. Doubleaveraging the Navier-Stokes equations produces the double-averaged momentum equations, whose projection on the streamwise x-direction is (Nikora et al. 2001): where d = − f ⟨ũ xũz ⟩ and t = − f ⟨u � x u � z ⟩ are called the dispersive (or form induced) and turbulent stresses, respectively. f p,x and f v,x are called the pressure drag and viscous drag on the solid element surfaces. v is the viscous stress.
Here, the turbulent stress t and the dispersive stress d have to be estimated from experiments. To that end, we must first measure the spatial disturbances ( ũ x ,ũ x ) and the fluctua- in a specific ROI. The velocity disturbance at any position can be esti- is the local time averaged velocity and ⟨u i ⟩ is the doubleaveraged velocity in a thin layer parallel to the mean bed surface at the mesoscopic scale. (The thin layer width is controlled by the vertical resolution of velocimetry measurements.) In the x-direction, we have U x = ⟨u x ⟩ , and if the flow is unidirectional at the mesoscopic scale, we also have ⟨u y ⟩ = ⟨u z ⟩ = 0 . A necessary condition for the flow to be considered two-dimensional is that these equations can be verified experimentally, and that the flow depth is uniform in the x-and y-directions. For further information on the fundamentals of double-averaging, the reader is referred to (Nikora et al. 2007).
In the surface layer (i.e. for z > z rc where (z) = 1 ), the double-averaged momentum equation is: In the permeable bed below the roughness crest (i.e. for (z) < 1 ), these assumptions are no longer valid because of drag interactions.

Constraints on laser sheet displacement when measuring using scanning
Section 2.3 presented the scanning methodology used to detect bead positions and acquire porosity profiles. Within this procedure, fluid velocity measurements can also be measured from the collected images. Coupling tomography and velocity measurements simultaneously are at the core of the PIV-RIMS method developed in this paper. This choice has the advantage of reducing experiment duration since it does not require performing several PIV measurements with fixed laser sheets. However, it requires imposing constraints on the laser sheet velocity. In the following, the sampling rate is denoted by f s ; thus, 1∕f s is the interval between two pairs of images from which PIV is performed. Δ t is the interframing time which is set to perform quality measurements of maximum velocities. In general, 1∕f s ≠ Δ t because multiple processing plans exist for a given set of images (images in a video may be skipped). A special case is 1∕f s = Δ t when displacements are calculated from consecutive images in a regularly spaced image sequence. Three constraints on the moving laser sheet speed velocity V MLS have been identified. The first regards the laser sheet thickness, and the second is related to the flow's spatial variations. When the flow is laminar, the first two conditions are sufficient to set V MLS and obtain the mean flow characteristics. However, when turbulence occurs, a third and generally stronger condition on V MLS must be considered to capture the local turbulent statistics. These three conditions might be helpful for managing the scanning procedure and optimise V MLS and f s for a given Δ t in various settings.

Laser sheet thickness constraint
In PIV, the laser sheet thickness w LST imposes that measurements are spatially averaged over a thin volume. With the translation in the transverse direction at a constant velocity V MLS , the displacement between two measurement slices is then V MLS Δ t . If this distance is too long, it will increase the number of uncorrelated particles. One first requirement is thus to ensure that V MLS Δ t is negligible compared to w LST . This gives the following constraint on V MLS : In the present configuration, Δ t = 1∕420 ∼ 2.3 ms with w LST ∼ 1 mm giving w LST ∕Δ t ∼ 0.42 m/s. Assuming that "negligible" means at least 40 times lower, we can reasonably consider that uncorrelated data due to the displacement of the laser sheet for the given laser sheet thickness is negligible if V MLS < 1 cm/s. In other words, as soon as this condition is met, we can assume that PIV measurements are similar to measurements that would have been made using fixed measurement planes.

Spatial variability constraint
The second constraint concerns spatial flow variability and the sampling rate f s . As mentioned above, the laser sheet thickness w LST imposes a spatial averaging. Another length scale L u can be selected by the experimenter if the flow's spatial variability is higher. It can be defined by the Taylor scale or by the typical length scale of the medium (a fraction of the grain diameter, for instance). To understand this constraint, we can make an analogy with a flatbed photograph scanner, which requires an adjustment of the carriage velocity for a given digital resolution. This second constraint is given by: For our set-up, measuring velocities between each frame ( 1∕f s = Δ t case) and setting L u to its lowest boundary w LST gives V MLS < 420 × 0.001 ∼ 0.42 m/s. We see here that the first constraint (6) on the laser sheet thickness is more V MLS < f s L u restrictive. As a consequence, it is possible to reduce f s if desired by skipping images (to reduce the amount of stored data for instance). With V MLS = 1 cm/s, we can set f s = 20 Hz and still meet the condition (7). In this situation, we record velocity every 0.5 mm. In practice, f s can be set at a higher rate to obtain additional measurements on overlapping volumes and then reduce inaccuracy.

Turbulence constraint
The third constraint is imposed by turbulence: the time during which the camera monitors a given flow slice must be sufficiently long to obtain accurate estimates of turbulent statistics, e.g. for estimating local turbulent stress. A wait during a specific time interval T u ′ is required to access the second-order statistics. This gives the following constraint on the moving laser sheet speed: It is worth mentioning that T u ′ also depends on f s . To record turbulent statistics, as fast as possible, it will be required to set f s at a high value to obtain samples of correlated data giving a fast convergence for turbulence statistics. However, T u ′ is physically bounded by the presence of low flow frequencies. In practice, and for classic high-speed camera while f s could be equal to 1∕Δ t , it can be set at a lower value than 1∕Δ t without modifying the time T u ′ required for convergence. In the following, we set f s at 210 Hz while keeping the inter-framing rate at 1∕Δ t = 420 Hz.
Because no preliminary information is available on the characteristic time T u ′ , experimentalists must estimate its order of magnitude. Section 3.5 tackles this issue. The next section outlines the averaging procedures.

Scanning and averaging procedure
Once the laser sheet has been shifted, time averaging is defined as: where is any local quantity of interest, such as the velocity u x (t, x, y l , z) or the instantaneous shear stress u � x u � z (t, x, y l , z) . Time T MA is the moving-average time, that is, the time window over which time averaging is done to obtain average local flow and turbulence statistics. Measurement is taken at time T m , and thus, the time window is centred on it. We denote the laser sheet's position by y l (t) : y l = V MLS t + y 0 , where y 0 is its position at t = 0.
Time-averages are then space-averaged over the y-axis. Averaging over the period T MA at time T m = Y m ∕V MLS implies that the averaging is done over the length D MA = V MLS T MA around Y m (see Fig. 5). The following condition is then obtained by matching the time-and space-averaging: Note that there is only one measurement at time t and at position y l during the translation. As the laser sheet has a finite thickness (about 1 mm), this thickness has to be included in the length D MA . When using PIV-RIMS techniques, time and space dependencies are intertwined, and it is therefore crucial to check the procedure's reliability with great care by comparing the RIMS measurements with those obtained using a fixed laser sheet at Y m over a long time period.

Flow characteristics and evaluation procedure
To assess scanning performance, we conducted two runs using the hydraulic characteristics detailed in Table 2. The bed arrangement was the same in both runs, but the runs differed as follows: • In the first run, velocities were obtained using PIV and a fixed laser sheet (FLS) positioned at Y m,FLS = 25 mm. • In the second run, velocities were obtained using the PIV-RIMS methodology and the flume was scanned using an moving laser sheet from Y m,0 = 2 mm to Y m,f = 40 mm.
First, we estimated the lag time T u ′ required to obtain accurate time-averaged quantities from the fixed laser sheet measurements. A 20-s period gave a robust estimate of the turbulence statistics near protuberances, making it possible to estimate T u ′ . The velocity of the moving laser sheet V MLS could then be deduced from the constraints imposed by Eq. (7) and Eq. (8). Scanning performance was evaluated by comparing the mean velocities and turbulent stresses obtained using the RIMS methodology and those obtained by the fixed laser sheet. Table 3 summarises the information. Figure 6 shows a snapshot of the vertical and horizontal velocity components, their time-averaged fields and the resulting space-and time-averaged velocity profiles with the laser sheet fixed at Y m = 25 cm. The light green areas give an idea of the spatial variability of the time-averaged velocities. The same is then done for the velocity fluctuations, disturbances and turbulent stresses. The time-averaged vertical velocity field u z helps to understand why spatial averaging is useful (see Fig. 6b2): u z exhibited large spatial variability. At this local scale, the flow was neither uniform nor unidirectional. It was only at the mesoscopic scale, after appropriate space-averaging, that the flow could be considered uniform and unidirectional. As shown by Fig. 6b3, the space-averaged vertical velocity profile ⟨u z ⟩ was close to zero across the entire depth, confirming flow uniformity.

Temporal and spatial averaging measurements using the fixed laser sheet
The time-averaged turbulence intensities along x and z also showed substantial spatial variability (see Fig. 6c2, d2).
| values were observed behind protuberances due to their generation of turbulent wakes. | | u ′ z | | remained more homogeneous, but its magnitude was higher, not only in the turbulent wakes but also against bead front faces on top of the bed. Inside the permeable bed, the turbulent activity was negligibly small. This observation must be tempered owing to the errors made when measuring small velocities using PIV. Indeed, we had Δ | | u � i | | small ± 2 mm/s . The highest velocities at the free surface were similarly subject to greater inaccuracy owing to the difficulties in measuring displacements at the surface (see Fig. 6c3, d3). If vertical turbulence intensities were assumed to be zero at the free surface, then the observed fluctuation might have resulted from the inaccuracy at this elevation as we had Δ | | u � i | | high ± 5 mm/s. We found that the turbulent stress t roughly matched the depth-averaged momentum flux, as expected from the momentum balance equation (5) when the dispersive and viscous stresses can be neglected inside the surface layer. The spatial disturbance fields ũ z and ũ x also exhibited large spatial variability (see Fig. 7f2, g2), but contrary to the turbulence intensities, their peak values were observed around the beads rather than in their wakes. Note that the spatial disturbance could be positive or negative, which was not the case for the turbulence intensities. As observed in Fig. 7 − (f2), there were small zones in front of and behind beads where the horizontal velocity components were lower. In these zones, vertical velocity was more often oriented upward. As a result, the dispersive stress d = ⟨− fũxũz ⟩ was more often positive than negative at the interface (see Fig. 7h3). For the FLS setting, measurement was taken in a single flow slice, and the dispersive stress was negative on top of the bed. This feature was not observed when averaging a larger domain using the PIV-RIMS procedure, as shown in Sect. 5.1.1. The dispersive stress affected both the front and rear of protuberances, i.e. where the velocity deficit was significant (see Fig. 7h2).

Turbulence statistics
We studied the turbulence around the protuberance formed by one of the borosilicate beads on top of the permeable bed, as shown in Fig. 8a, by using PIV and the fixed laser sheet (FLS) aimed at eight points around that protuberance. The spatial variability around this protuberance can also be observed in the two-dimensional time-averaged statistics of Fig. 7.
For the measurement points located above the roughness crest (A1, B1, C1), turbulence was spatially homogeneous and of weak intensity. For the measurement points on the roughness crest (A2, B2, C2), the intensity of turbulence was higher and differences between measurements point were visible. Finally, for the lowest level, in the rough layer (A3, C3), there was large spatial variability in turbulence. The average velocity at point C3 was close to 0 and the signal-tonoise ratio was therefore very low, whereas for A3, upstream of the bead, the velocities were higher, with high turbulence intensity relative to the mean velocity. Because Mignot et al. (2009a, b) have thoroughly explored these flow types, we will not go into this topic further. Here, statistical analysis of the turbulence identifies the region where T u ′ was the largest. Figure 8b shows how the empirical error in the averaged velocity computation depended on the measurement's duration T. For the flow zones around point C3, we found that T u ′ had to be as long as 2 s to obtain a relative error lower than 10%. This was the strongest constraint to our continuous scan methodology.
In the configuration tested, and with T u � ∼ 2 s, the more restrictive condition was given by inequality (8) since Δ t ≤ 1∕f s ≪ T u � = 2 s . The maximum velocity required by the MLS to obtain satisfactory continuous scan measurements could thus be estimated using Eq. (8) at V MLS,max ∼ L u T u � ∼ 2 mm/s if L u was approximated by d p ∕2 ∼ 4 mm. Ideally, it would be better to set the length scale value at its lower bound given by the laser sheet thickness w LST . However, it would impose a long scanning procedure incompatible with our experimental conditions. As we will see in the next section, V MLS = 2 mm/s was satisfactory to determine flow statistics.

Evaluation results
After the scanning velocity was determined using Equation 8, the MLS run was performed by setting V MLS = 2 mm/s. Figure 9a1, c1 shows the time-averaged velocity and turbulent stress fields at position Y m = 25 mm from the MLS. Figure 9a2, c2 compares the resulting profiles (averaged along the x-direction) using the FLS and MLS procedures, respectively. Figure 9b1, d1 shows the absolute differences observed between the FLS and MLS procedures carried out on the field and on the averaged profiles. Averaged velocity and turbulent stress at laser sheet position Y m showed good matches between the FLS and MLS procedures.
Error estimates were obtained by subtracting the FLS velocity and turbulent stress profiles from those acquired during the MLS experiment at Y = 25 mm (see Fig. 9b2, d2).
Various times T MA were tested and, as observed in Fig. 9e1, e2, the smallest errors were obtained for T MA = 2 s, which corresponded to the prediction made in Sect. 3.5.3. If T MA was shorter or longer than 2 s, the error increased.

The sidewall's influence on the flow
To demonstrate the potential of the RIMS method, we present an investigation of the sidewall's influence on the flow. Figure 10 shows a three-dimensional reconstruction of the flow, that is, the horizontal velocity component on the wall of a cuboid that represents the ROI. The fluid's velocity u x increased with y as measurements were taken further away from the sidewall. This increase can also be observed in Fig. 11, where u x has been averaged in the x-direction and plotted for different z � =0.8 . We found that the flow region which felt the sidewall's influence least was Y = 10 mm away from the wall. This demonstrates that measurements taken at less than 5 mm from the wall were strongly affected by it. These average velocity estimated using n = 700 samples of duration T against the empirical average calculated over T tot = 20 s. The uncertainty is under 10% after 2 s Fig. 9 Comparison of estimated profiles at position Y m = 25 mm using the fixed and moving laser sheet methods. The error is thus defined by where is the velocity profile U x or the turbulent stress t . T MA is the averaged time or, equivalently, the distance D MA framing the position Y m . These results show that error is low for both velocity measurements and turbulence statistics when T MA ∼ 2s observations provide a posteriori grounds for using the index matching method, exploring the flow at a sufficient distance from the sidewall and obtaining profiles that can be used to evaluate the different contributions to the double-averaged momentum equation.
Interestingly, in Fig. 11, the velocities close to the free surface (z � =0.8 = 9 mm) and the wall (Y < 7 mm) were lower than velocities in deeper positions ( z � =0.8 = 2 mm and z � =0.8 = 4 mm4 ) but at the same distance from the wall. This phenomenon can be understood by noting that shear and dispersive turbulence were stronger next to the bed. The resulting mixing processes actively convected momentum from the middle of the flume to the sidewall at these depths, whereas near the free surface, the momentum transfer was of lower magnitude.

The bed arrangement's influence on reproducibility
When using the PIV-RIMS methodology under similar flow conditions (slope and flow rate) and averaging measurements between Y = 10 mm and Y = 40 mm (to avoid sidewall influence), consecutive double-averaged velocity profiles were similar, thus indicating the experimental procedure's good reproducibility as long as the bed was not rearranged (see Fig. 12).
As the mesoscopic scales ( ∼ 8 cm in the x-direction and ∼ 3 cm in the y-direction) used in the double-averaging were not much larger than the bead size (d p ∼ 1 cm) , our measurements suffered from finite size effects. In other words, we could not guarantee that porosity and averaged velocity  7). b Frontal view of the flow, sliced along y. This view enables us to appreciate the sidewall's influence. Y is the distance from the sidewall. The different colour of the beads (purple vs pink) represents the two different diameters (7 and 9 mm). Z and X are arbitrarily referenced Fig. 11 Sidewall influence. The horizontal velocity profile has been averaged along the streamwise x-direction and plotted for different elevations z � =0.8 = z − z =0.8 as a function of their distance from the sidewall Y Fig. 12 Comparison between two consecutive runs using the PIV-RIMS procedure, with each run using identical bed structure and flow characteristics [left] Velocity profiles [right] turbulent stress profiles. Velocity profiles look very similar while a slight difference is observable in turbulent stress profiles profiles were insensitive to changes in the bed arrangement. We now take a closer look at this issue. Figure 13 compares ten porosity and velocity profiles. These profiles were measured using a constant flow discharge (q f = 0.30 dm 2 ∕s) and varying slopes ( 1% , 2% and 4% ). Initially, the bed was randomly mixed and flattened using a ruler. The flow depth was still h s = 5 cm. The ROI was located at a distance g = 90 cm from the outlet (see Fig. 14). How flow uniformity depended on this value is addressed in the next section. As we had found that slight variations in the slope could affect these profiles, we reset the slope's incline before each run and measured it to within 0.1%.
To compare vertical velocity and porosity profiles, we first needed to fix the origin of the vertical axis. There is no standard procedure for doing this with rough beds. We addressed various possibilities. The roughness crest was unsuitable because it created significant scatter between profiles: at the mesoscopic scale, z =0.99 was highly influenced by individual grains that were slightly higher than the average bed level. The origin had to be fixed at a bed height where the scatter between the porosity profiles was minimal. We found The modified coordinate is givenz =0.8 (see Fig. 13), that is, at 0.3d p below the roughness crest, which was the shift that Voermans et al. (2017) obtained using the porosity inflection method, and which was close to what Nezu and Nakagawa (1993) prescribed. This similarity with the findings described by Voermans et al. (2017) made it possible to compare their results and ours. The Reynolds numbers based on this choice were slightly different to those estimated from the height of the roughness crest, as the flow depth h f was computed from z =0.8 , which is above z =0.99 .
In the roughness and surface layers, the porosity profiles plotted in Fig. 13 showed a similar pattern from one experiment to another. Slight differences were, however, observed near the roughness crest. Profiles in the subsurface layer (located at z � =0.8 < −0.5d p in all runs) showed more scatter. This was the consequence of the ROI's finite size. The porosity profile tended to the packed bed porosity b ∼ 0.4.

Fig. 13
Evaluation of reproducibility with different bed structures. a Velocity profiles (continuous and dash-dotted lines) and porosity profiles (dotted lines) for different slopes but using a constant flow discharge q f = 0.30 ± 0.015 m 2 ∕s . The modified coordinate is given by z � =0.8 = z − z =0.8 . b A zoom in on the roughness and subsurface velocities Fig. 14 Permeable outlet condition to ensure a subsurface flow. Measurements must be taken along a sufficiently long distance g to ensure that this condition's boundary effect is negligible Although data exhibited a significant scatter, spatial averaging scales were sufficiently large to discriminate slope influence and bed arrangement effect. (We observed a low signal-to-noise ratio.) Based on this observation, we expected that the domain could we viewed as a representative elementary volume (REV). A larger domain would have improved similarities between runs with rearranged beds.

Uniformity: the influence of the permeable grid
At the flume's upstream end, honeycomb-shaped straighteners stabilised the inflow created by the constant head tank. Downstream of these straighteners, the flow ran over the bi-dispersed borosilicate beads (as shown by Fig. 3a). At the flume outlet, a permeable grid located at x g let the flow seep out of the bed in such a way that the subsurface flow occupied the bed's entire height (see the enlarged view of the flume in Fig. 14). Most experimental investigations of supercritical flows neglect the downstream boundary's influence, but in our experiments involving high bed-permeability, we observed that the downstream boundary condition affected a long section of the flume's length. We thus believe that it is essential to take this influence into account. Indeed, when, for instance, the permeable grid created excessive runoff from the granular bed, the surface flow got into the bed upstream of the grid, causing a decrease in the flow depth over a certain length of the flume g from the grid. The flow depth was then non-uniform over a more or less long part of the flume.

The permeable grid's influence on the subsurface flow
Within the Darcian framework (Darcy 1856), a quantitative estimation of g , i.e. the distance from the grid, where the ratio between the estimated surface flow rate and the theoretical surface flow rate in a uniform situation is , is estimated and given in "Appendix 2". It results in the following relation : With BAE = 3 × 10 −6 m 2 ∕s and d p = 8 mm , permeability was estimated using the Kozeny-Carman equation K = 3 d 2 p 180(1− ) 2 ∼ 6.32 × 10 −8 m 2 . The granular bed depth was set at h s = 0.05 m . With i = 2% (the experiment's average slope) and = 0.8 , we obtained a distance 0.8 g = 0.68 m . This method predicted that the permeable grid's domain of influence was fairly long relative to the flume length (∼ 2 m) . To control both the subsurface and surface outlet flow rate, we added a buffer medium (BM) with a permeability higher than that of the granular bed at the flume outlet. Moreover, the above prediction is made within the Darcian framework by considering only a viscous drag in the bed. In reality, we suspect form drag to also play an active role, especially in the roughness layer. Taking this additional drag into account can lead to reduced distance predictions. In the next subsection, we will take a closer look at the occurrence of a uniform flow along the flume.

Flow uniformity
To verify that the flow was uniform for x ≤ x g − g = 90 cm , we conducted experiments by varying the g length from the outlet. The longest distance was g = 110 cm , and for this value every position along the flume length was within the boundary's domain of influence. The shortest distance was g = 60 cm . Figure 15 shows that the outlet condition affected the flow for g distances as long as 60 cm for i = 1% . Higher velocities were measured in both the surface and subsurface flows, whereas the depth was lower than the averaged velocity profiles at g = 90 cm . This was expected from Eq. 21, where any increase in the subsurface flow rate caused the flow depth to decrease. For i = 4% , differences between profiles could not be statistically attributed to the outlet condition's influence, given the uncertainties and Fig. 15 Velocity profiles for various g values with which to evaluate the flow uniformity condition along the flume. The dashed-dotted lines represent the averaged velocity profiles at g = 90 cm. The error bars show the deviations from the averaged profiles due to the modification of the bed structure. They represent the 95% confidence interval. U x (z) is the standard deviation at z calculated from the profiles shown in Fig. 13. The continuous lines are the profiles measured at g = 60 cm, and the dashed lines were measured at g = 90 cm. The bottom-right inset plots the same profiles using a logarithmic scale to emphasise the marked differences at low velocities noise induced by the bed arrangement. This analysis suggests that a nearly uniform flow was reached at g = 90 cm because the differences between the profiles at g = 110 cm and g = 90 cm were not statistically significant.

Preliminary observations based on the PIV-RIMS method
As stated above, the main advantage of the PIV-RIMS methodology is that it averages flow quantities over a thin layer parallel to the bed, as prescribed in the double-averaging framework. Figure 13 shows double-averaged velocity profiles from various different bed arrangements and flume inclination.

Slope and averaged velocities
In the velocity profiles shown in Figs. 13 and 15, the slope was increased from 1 to 4%. As the slope was increased, we found that fluid velocities increased in all flow layers, but the U x increment varied differently according to the layer considered. As the slope was increased from 1 to 4% (see inset of Fig. 13), we found that the averaged subsurface layer velocities were multiplied approximately fourfold, whereas free surface velocities were multiplied by 1.5. This difference reflected the various mechanisms at play in those layers: flow through the porous medium was controlled by drag forces, whereas surface flow was mostly driven by the vertical momentum transferred by turbulence. Moreover, (Breugem et al. 2006) suggested that all flows involving a sediment-fluid interface exhibit an inflection point in their velocity profiles. This suggestion was confirmed here and is also consistent with the observations by Voermans et al. (2017). Figure 16 shows the dispersive and turbulent stresses obtained using the FLS and PIV-RIMS methodology. One significant difference between the two experimental procedures-use of FLS or MLS-was observable in the measurement of dispersive stresses. The turbulent stresses computed by both procedures showed similar behaviours. Using the FLS, the averaging procedure was done for a single slice at Y m = 25 mm, eliminating the possibility of collecting information at other y positions. Using PIV-RIMS, the laser sheet moved along y and the profiles were averaged over x and y, thus reflecting flow variability in the transverse direction. As observed in previous studies (Voermans et al. 2017;Fang et al. 2018), the dispersive stress exhibited a positive trend at the interface, with maximum dispersive stress located just below the roughness crest (here, at z � =0.8 ∼ 0 ).

Dispersive and turbulent stresses
The turbulent stress maximum was located slightly above the roughness crest, at z � =0.8 ∼ 0.3d , and rapidly decreased with decreasing z ′ (or, equivalently, with increasing porosity).
Turbulent and dispersive stresses rapidly dampened in the subsurface layer, i.e. for depths below z < z t , where bed porosity reached ( z � =0.8 < −0.5d in Fig. 16). Flows in the deepest layers were indeed essentially controlled by drag forces on grains. The roughness layer, that is, the transition zone, where porosity varied sharply from z rc to z t ( 0.3 d > z � =0.8 > −0.5 d in Fig. 16), was the zone which presented substantial vertical momentum exchanges. These exchanges resulted from either turbulence or dispersive effects.

Concluding remarks
The present article presented a PIV-RIMS technique for measuring averaged flow variables at the mesoscopic scale (velocities, stresses and porosity) as part of the doubleaveraging approach. The technique was applied to a turbulent unidirectional flow over a porous, coarse-grained bed. Combining laser scanning and iso-index techniques (RIMS) made it possible to obtain accurate porosity profiles (z) . Coupled with the PIV processing, this technique enabled us to determine the velocities in the surface, roughness and subsurface layers of the fluid.
The PIV-RIMS methodology is a fast and convenient experimental procedure, but it requires adjustments to the Fig. 16 Comparison of the dispersive and turbulent stresses obtained using either a fixed laser sheet (FLS) or the PIV-RIMS methodology. In contrast to using the FLS procedure, PIV-RIMS captures the variability of interactions in the transverse direction y. The resulting averaged profiles provide a better representation of the profiles at the mesoscopic scale experimental parameters: the velocity of the MLS V MLS has to be slow enough to extract the flow's spatial and temporal variability. To measure mean flow properties far from the flume sidewalls influence, we computed velocity profiles by averaging the flow between Y = 10 mm and Y = 40 mm.
As the dimensions of region of interest were constrained by the flume, measurements were sensitive to bed arrangement between runs. However, reproducibility tests were conducted successfully, allowing for measurements of the slope influence on subsurface flows. We also found that we had to place the region of interest sufficiently far from the flume outlet-at g = 90 cm-to ensure flow uniformity on average.
Preliminary observations revealed the roughness layer's crucial role in the transfer of vertical momentum. These observations were thus in contrast with the assumption commonly used in most extant models-that there is a discontinuous porosity profile at the bed-flow interface. A valuable alternative might be to consider the roughness layer and model the continuous variations in the velocity and porosity profiles. In a near future, the measurement system could be rotated to see the flow from below and probe transverse and streamwise velocities. The system could also be improved by combining the stereoscopic PIV technique to obtain in one scanning procedure the three velocity components in a three-dimensional domain (3D3C-PIV).
Our results open up another avenue for refining closures involved in current models working at the mesoscopic scale. These data have been used in Rousseau (2019) PhD thesis where various modelling assumptions for low relative submergence flows were tested with an application on mountain river processes. In general, these data should contribute to a better understanding of flow resistance or transport processes in various settings (e.g. mountain rivers processes, hyporheic processes or heat transfers). These results could also be useful to test and validate large-eddy-simulation models (Fang et al. 2018;Lian et al. 2019) or models based on double-averaged or Reynolds averaged momentum equations (as in Maurin et al. (2018)

for instance).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.

Appendix 1: Image velocimetry processing
This Appendix details the three principal steps of the image velocimetry algorithm used to yield the velocity field from consecutive images. The test on PIV challenge Case A (Kähler et al. 2016), as well as the internet address needed to access to the algorithm's detail, are provided at the end of this Appendix.

Pre-processing
For a given laser sheet position Y m , a mask is generated from the bead positions to restrict measurement to the interstitial flow zones and the surface flow. Similarly, the fluctuating fluid/air interface is detected in order to mask the upper part of the frame (see Fig. 17a).
The Contrast-Limited Adaptive Histogram Equalisation (CLAHE) algorithm (grid size = 32× 16 px, clip limit = 8) is used to improve image contrast. Hot pixels (with constant high-intensity values) may be present during the recording, as may local and temporary (long duration in comparison with the particle displacement) hot spots due to reflection. To solve this problem, a background removal procedure is performed by subtracting the average frame.
Before any velocity measurements, the Good Feature to Track (GFT) algorithm selects features that maintain good contrast, i.e. that are able to provide accurate velocity estimates (Shi 1994) (see Fig. 17b). This pre-selection has two advantages: first, by discarding low-quality points it diminishes errors, and second, it decreases the number of potential interrogation windows, thereby making the algorithm more efficient. Using classic PIV methods, the entire domain is usually computed within regularly spaced interrogation windows and low-quality measurements are generally discarded using post-processing methods. The present study's method avoids processing those zones with a low signal-to-noise ratio. After this step, the points i = (x, z) T i are selected and velocimetry processing is launched.

Velocimetry processing
The velocimetry algorithm measures the local (regionbased) optical flow by means of a pyramidal implementation of the Lukas-Kanade method (Bouguet 2001) (see Fig. 17c). This method minimises the square of the Displaced Frame Difference. To better understand the equations involved in this algorithm and its link to classic PIV, the procedure is detailed below, based on the papers by Heitz et al. (2010) and Liu and Shen (2008).
Given a position = (x, z) T in the image and the intensity function I( , t) of the image field, the velocity field is denoted as The Optical Flow Constraint (OFC) equation representing the brightness constancy can be written as Equation 13 is the linear formulation of the matching formula between two consecutive images and is also known as the Displaced Frame Difference: where ( ) denotes the displacement field between two images. With the Lukas-Kanade method, the displacement field between two consecutive images is determined by minimising the square of the Displaced Frame Difference model where W( ) is the interrogation window around the point of interest. Since I( , t) is independent of , Equation 15 is equivalent to: Equation 16 shows that the minimisation of the square of the Displaced Frame Difference includes the correlation between two consecutive images. The displacement field estimated using this method is thus equivalent to the displacement field obtained using classic PIV if the quantity ∑ ∈ ( ) I( + ( ), t + t) 2 does not depend on . Classic cross-correlation techniques implicitly use this assumption, but it is locally strengthened when small interrogation windows or large velocity gradients are considered. This is probably why this method works well for the current study problem, where small pore sizes limit the windowing.
The pyramidal application of the algorithm aims to increase its dynamic range, i.e. deal with significant pixel motion. The pyramid refers to the successive low-pass filtering and sub-sampling of the image sequence. The levels of the pyramid (1 2, 3, ...) represent the number of passes and the resolution of the image for the first pass on which the Lukas-Kanade velocimetry method is executed. For example, if the image has a resolution of 400×400 px and the pyramid has two levels, the first image has a resolution of 100×100 px. In this image, the pixel motion is smaller, and the Lukas-Kanade method (with the same window size) measures the overall movement to introduce a shift for the second pass. This methodology is the equivalent of the iterative multi-grid method commonly used in fluid mechanics (Scarano and Riethmuller 1999). For this experiment, the pyramidal Lukas-Kanade method is parametrised with a 16 × 16 px window and three pyramidal levels.
At the end of this step, the velocity is obtained for each of the selected points i = (u x , u z ) T i (see Fig. 17c).

Post-processing and interpolation scheme
The final step consists of an interpolation process to obtain a velocity field from the isolated points where velocity was known, filling the gaps where the image quality was poor or where the number of seeding particles was too small. This step is commonly performed when using particle tracking velocimetry (PTV) algorithms, but it is computationally expensive. Recent improvements in the Visualisation ToolKit (VTK) library allow use of a tree-like data structure to partition the 2D space and create buckets (methods that are commonly used in 3D graphics or 3D game engines). The search for the points or closer neighbours is then more efficient. Within the Darcy framework, the subsurface layer flow is not expected to exhibit linear behaviour at the outlet, where velocities increase. However, the Ergün equation's quadratic term usually decreases the permeability, and the linear approximation has the effect of overestimating the flow inside the porous bed at the outlet. Thus, the increase in the total subsurface flow discharge q f ,h SSL is given as a function of g : At this point, it is observed that as g → +∞ , flow discharge in the bed tends to its expected steady value q 1 f ,h s = Kg h s i and the steady surface flow is given by q 1 f ,SL+RL = q f − q 1 f ,SSL . Equation 19 involves h f ( ) which is non-uniform along x. To resolve this equation, a relation between surface level and subsurface flow discharge must be provided. It is quite complex due to the non-uniformity of both surface velocity and depth. Instead, h SL is supposed to be negligible with respect to h SSL + h RL . This assumption seems reasonable since h f is about 0 at the outlet condition. Also, we must scale h RL and h SSL . As observed experimentally, h RL ∼ d p and the subsurface layer thickness are given by h SSL ∼ h s − d p , where h s is the initial total sediment depth fixed manually. The next step produces the distance g , where the surface flow decrease is negligible. With the condition that q f ,SL > q 1 f ,SL , where is the quality coefficient (that should be close to one to obtain a nearly uniform flow), we obtain the following equation: The height , above which this condition is verified, is thus provided by: