Momentum balance in flows over mobile granular beds: application of double-averaging methodology to DNS data

ABSTRACT This paper presents the application of the double-averaging methodology to actual data obtained for the mobile-bed conditions from direct numerical simulations using an immersed boundary method. The dimensions of the computational domain resemble those for open-channel flows with small relative submergence. The domain bottom is a plane covered with one layer of hexagonally packed, single-size spheres fixed to the bed. The fixed particles are covered by 2000 mobile particles of the same size that are free to move. Two simulation scenarios at distinctly different values of the Shields parameter are studied representing a fully mobilized and partially mobilized granular bed. The effects of the averaging time and the averaging domain size and shape on the double-averaged flow quantities are identified first. Then, the data analysis focuses on the detailed assessment of the key terms of the double-averaged momentum balance equation formulated for mobile-bed conditions. The paper demonstrates that the double-averaging methodology provides an efficient reduction of the massive datasets produced by the fully resolved simulations to a manageable number of physically meaningful double-averaged quantities, which should help devising closure strategies for modelling mobile-bed flows.


Introduction
Turbulent flows over mobile beds are ubiquitous features of many environmental and engineering systems. They largely determine systems performance in terms of mass, momentum, and energy transport and therefore have attracted significant attention, particularly in environmental applications (Yalin & Ferreira da Silva, 2001). However, most applications and modelling approaches for predicting bed stability and sediment transport rates, starting from the seminal papers on the topic of Meyer-Peter and Müller (1948) and Bagnold (1956Bagnold ( , 1966 remain to be grounded on the 80-year old framework of Shields (1936), as reviewed by Buffington and Montgomery (1997). The approach of Shields (1936) is based on a simplified balance of forces acting on bed particles considering the ratio of the destabilizing forces (drag, lift) to resistive forces (gravity, friction), today known as the Shields parameter. The critical value of this parameter, which separates particle stability from entrainment by a turbulent flow, was found to depend on the particle Reynolds number. The Shields parameter is typically used for the assessment of both bed stability and bed-load transport, which are assumed to be a function of the excessive bed shear stress, i.e. the difference between its actual and critical values. Although this concept is known to be problematic (Seminara, Solari, & Parker, 2002), the approach is still widely used in hydraulic applications (Garcia, 2008) in spite of large errors in estimates of bed stability and bed-load transport, often exceeding 100%. This high level of uncertainty motivates new studies into mechanics of bed-load transport and its effects on the overall flow dynamics (e.g. Bathurst, 2007). The list of key issues waiting for clarification is long and probably should start with rigorous analytical frameworks for sediment motion and flow dynamics over mobile beds, which are currently under intensive development (e.g. Ancey & Heyman, 2014;Nikora, Ballio, Coleman, & Pokrajac, 2013). These frameworks should provide a sound platform for studying the multifaceted mechanisms of sediment entrainment, transport, and deposition. An important part of this development relates to advancing the current knowledge of grain-scale processes so that they can be rigorously accounted for by large-scale models used in applications. Linking sediment mechanics with flow dynamics through coupled hydrodynamic equations rather than employing empirical correlations or intuition seems to be the most appropriate strategy to address this issue.
The development of the ideas on grain-scale processes in sediment transport research remains slow due to the lack of high resolution data in relation to both time and space domains. Although experimental technologies for obtaining such datasets are emerging (e.g. Campagnol et al., 2013), it will take some time before they are used in full for revealing physical mechanisms and testing theoretical predictions. Direct numerical simulations (DNS) with spatially resolved solid-liquid interfaces using the immersed boundary method (IBM) as proposed by Kajishima, Takiguchi, Hamasaki, and Miyake (2001), Uhlmann (2005), or Kempe and Fröhlich (2012a) may serve as an alternative or complimentary approach, providing a huge wealth of data (Chan-Braun, García-Villalba, & Uhlmann, 2011;Kidanemariam, Chan-Braun, Doychev, & Uhlmann, 2013;Shao, Wu, & Yu, 2012;Vowinckel, Jain, Kempe, & Fröhlich, 2017). This approach, however, requires an appropriate representation of the physical configuration targeted. One important issue in this respect is the use of a suitable particle collision model to represent particle-particle interactions for dense suspensions, i.e. for collision-dominated and contact-dominated particle transport (Kempe & Fröhlich, 2012a;. Another issue is the identification of the domain size to represent the typical length scales of the flow and the particle structures Vowinckel, Kempe, Fröhlich, & Nikora, 2012).
Large computational domains result in very costly simulations and an enormous amount of data that need to be convoluted into a practical set of parameters. In other words, effects of grain-scale processes need to be up-scaled to a much larger scale of engineering applications dealing with a wide variety of sediment transport problems. This up-scaling can be based on the time-and space-averaged hydrodynamic equations applicable for both flow regions within and above the fixed and/or mobile granular beds following an approach known as the double-averaging methodology (DAM) (Nikora, Ballio, Coleman, & Pokrajac, 2013;Nikora, McEwan, et al., 2007;Nikora, McLean, et al., 2007). The DAM stems from the terrestrial canopy aerodynamics (e.g. Finnigan, Shaw, & Patton, 2009;Raupach & Shaw, 1982). In recent years, it has also been successfully employed for studies of open-channel flows over fixed rough beds, examples of which are given in Nikora and Rowiński (2008), Franca, Ferreira, and Lemmin (2008), Mignot, Barthelemy, and Hurther (2009) , Ferreira, Ferreira, Ricardo, and Franca (2010), Yuan and Piomelli (2014). Among other findings, these studies highlight the potential importance of additional (dispersive) fluid stresses, which appear due to roughness effects that introduce local heterogeneities in the time-averaged flow.
The goal of the present paper is to compute, for the first time, the double-averaged momentum balance of the fluid phase in flows within and above mobile granular beds composed of monodisperse, cohesionless particles of a finite size using highly resolved DNS data. This information will help to identify the most important terms and their relevance when developing models for mobile-bed conditions. The DNS data were obtained in previous work of the authors  for two situations: one well above and another well below the particle stability threshold for incipient motion characterized by the critical Shields parameter.
The paper is structured as follows. First, the numerical method and the computational set-up used in this study are outlined, followed by a description of the simulation scenarios. Then, a summary of the DAM is provided, highlighting the essential features of the method as applied in the present work. In this approach, temporal and spatial averaging domains are key factors and therefore their selection is discussed next. The main body of the paper includes a detailed assessment of the momentum balance for the studied scenarios applying the double-averaging framework. This results in new findings highlighting the significance of particular terms of the momentum balance relevant to the closure problem in the development of up-scaled numerical models for mobile-bed flows. A complementary paper by Vowinckel, Nikora, Kempe, and Fröhlich (submitted to JHR, 2016) provides additional information on spatially averaged momentum fluxes determined on the basis of the method presented here.

Immersed boundary method
In the case when particles are larger than the Kolmogorov scale, point-particle approaches cannot be used without additional modelling, which involves well-known uncertainties. Such point-particle approaches, for example, particularly suffer from the uncertainties in the required empirical correlations for drag and lift, which become large if many particles are close together or if particles are colliding (Soldati & Marchioli, 2012). Bed-load transport is characterized by high volume fractions of the solid phase and very dense particle clusters close to the sediment bed. This raises the need for fully-resolved simulations with a four-way coupling of the flow and the disperse phase (Balachandar & Eaton, 2010). In the present study, appropriate four-way coupling is realized using the enhanced IBM developed by Kempe and Fröhlich (2012a), which extends the method of Uhlmann (2005).
The DNS approach, employed in the present study, allows the simulation of a large number of fully-resolved mobile particles in a viscous fluid. It is fully described in Kempe and Fröhlich (2012a), so only a brief description is given here. The continuous phase is governed by the unsteady three-dimensional Navier-Stokes equations for incompressible fluids: where τ is the fluid stress tensor, u = (u 1 , u 2 , u 3 ) T = (u, v, w) T is the velocity vector with velocity components u, v, and w in streamwise, wall-normal, and spanwise directions, respectively, ρ f is fluid density, p is pressure, ν f is fluid kinematic viscosity, f IBM is a source term introduced by the IBM as detailed below, and I is the identity matrix. The term f represents an artificial volume force driving the mean flow; it is adjusted in each time step to enforce a given bulk Reynolds number. Note that in many typical hydraulic situations the momentum is supplied by gravity g = (g sin α, g cos α, 0) T , where g is the gravity acceleration and α is the bed slope angle. To give an example, the bed slope is of the order of α ≈ 0.01 for a similar set-up as in this study (Cameron, Nikora, & Coleman, 2008). In Eq. (1) the gravity effect is replaced with the volume force f = (f x , 0, 0) T . The wallnormal component f y is not accounted for in the fluid solver because the hydrostatic pressure contribution is not relevant to the problem under consideration. This is standard practice followed in all numerical studies of this type of configuration (e.g. Kidanemariam et al., 2013;Papista, Dimitrakis, & Yiantsios, 2011;Shao et al., 2012). The spatial discretization of Eqs (1) and (2) is performed by a second-order finite-volume scheme on a fixed Cartesian grid with staggered positions of the variables. Pressure is stored at the cell centre of the basic grid and velocity components in the centres of those cell faces to which they are orthogonal. The convective terms are treated with an explicit Runge-Kutta three-step method and the viscous terms with an implicit Crank-Nicolson scheme. The divergence-free velocity field at the end of each Runge-Kutta sub-step is obtained by the solution of a Poisson equation for the pressure correction and subsequent projection. A full description of the numerical scheme can be found in Kempe and Fröhlich (2012a) and Kempe et al. (2014).
Discrete marker points are used to represent the fluid-solid interface as illustrated in Fig. 1. The method introduces a source term f IBM (right-hand side of Eq. (1)) in the vicinity of the interphase boundaries to connect the motion of the particles with the Figure 1 Particle surface with 315 discrete marker points and three slices of the Cartesian Eulerian background grid representing the fluid covering the entire computational domain liquid phase. The transfer of quantities, like force and velocity, between Eulerian points and Lagrangian points, i.e. between fluid points of the regular background grid and points on the surface of the particle, is performed by interpolation and spreading operations via a weighted sum of regularized Dirac delta functions. The source term f IBM is computed in such a way that the no-slip condition at the particle surface is fulfilled. The fluid velocity is set equal to the particle surface velocity at each individual marker point on the interface employing the direct forcing approach of Mohd-Yusof (1997) and internal loops over the forcing (Kempe & Fröhlich, 2012a). A major advantage of such an IBM is the avoidance of permanent re-gridding and mesh adaptation in the case of mobile particles, as would be required for approaches with body-fitted meshes. Hence, the method is very well suited for the simulation of a large number of mobile particles .
The motion of spherical particles is described by ordinary differential equations for the translational velocity u p of a particle: and for its angular velocity ω p : where m p is the particle mass, ρ p is the particle density, n is the outward-pointing normal vector at the surface p of a particle, g represents gravitational acceleration, r = x − x p is the position vector of a surface point with respect to the centre of mass of a particle, and I p = 8πρ p R 5 p /15 is the moment of inertia of a spherical particle of radius R p . Particle-particle interaction is accounted for by source terms on the right-hand sides of Eqs (4) and (5), where F p represents the collision forces and M p is the moment due to particle interactions addressed in the following sub-section. The three-step Runge-Kutta scheme used for the fluid is also employed for the time integration of the governing equations of the particle motion, and the same time step is employed for both. As mentioned above, the gravity effect is not directly presented in the fluid solver (Eq. (1)), but is explicitly accounted for in the particle motion in Eq. (4) as g y = |g| appears in the buoyancy term.

Particle collision modelling
Typically, bed-load transport involves particle-particle collisions and contacts engaging both moving and immobile particles. As the particle collisions are key mechanisms of the momentum exchange, their adequate representation is of major importance in simulations of mobile-bed flows. The present simulation employs the adaptive collision model (ACM) recently proposed by Kempe and Fröhlich (2012b) and extended for multiple simultaneous collisions in Kempe et al. (2014). It accounts for unresolved lubrication forces, normal forces during direct contact of surfaces, and tangential contact forces, without affecting computational performance. The central idea of the ACM Table 1 Physical parameters common for both simulation scenarios, where particle Reynolds number D + = u τ D/ν f , bulk Reynolds number R b = U b H /ν f , and friction Reynolds number R τ = u τ H /ν f represent reference values obtained for the clear-water flow (with no moving particles); u τ is the shear velocity, U b is the bulk flow velocity, D the particle diameter, ν f the fluid kinematic viscosity, H the flow depth above the crests of fixed bed particles, N p ,fix the number of fixed bed particles, and N p ,mob is the number of mobile particles is to stretch the duration of contact in time, thus avoiding excessive time step reduction when a collision occurs in the simulation, while still maintaining physical realism by a judicious optimization process. The cited references contain a detailed description of the model as well as numerous validation studies for particle-particle and particle-wall collisions. The ACM provides a substantially improved representation of collision and contact compared to other models, such as that of Glowinski, Pan, Hesla, Joseph, and Periaux (2001). This latter model is based on a repulsive potential and is usually employed with imposing a minimal distance between particle surfaces, so that mobile sediment particles have a somewhat higher protrusion into the flow and interact with the bed underneath in a different, less realistic way .

Computational set-up
A turbulent open-channel flow is considered in a rectangular computational domain with periodic boundary conditions in streamwise and spanwise directions, a free-slip condition at the top, and a no-slip condition at the bottom and the particle surfaces (Table 1, Fig. 2a). The sediment bed is composed of a single layer of fixed spheres of diameter D, in hexagonal packing, in such a way that the line of closest connections between adjacent particles is perpendicular to the mean flow. The origin of the wall-normal coordinate y is set at the crest of the fixed spheres, with the flow depth H defined as the distance from the particle crests to the domain top. The computational domain is = [0, 12 H ] × [−D, H ] × [0, 6 H ] with extents L x , L y , and L z in streamwise, wall-normal, and spanwise direction, respectively. The chosen shape and extents of the computational domain reflect the typical wavelengths of the coherent structures observed in turbulent flows over rough walls, which are l x ≈ 6H and l z ≈ 2H, respectively (Marusic et al., 2010;. Indeed, in  it was reported, based on a two-point correlation analysis of the present Yellow: fixed bed particles; white: mobile particles with |u p | < u τ ; black: mobile particles with |u p | ≥ u τ data, that the selected spanwise extent of the simulation domain serves well to display the relevant spanwise length scales. For the streamwise direction,  showed that the correlations do not decay to zero within the selected longitudinal extent of the domain. This is attributed to the effects of the fluid-particle interactions and particle clusters on the bed. As will be discussed in Section 3.1, heavy particles form streamwise-oriented ridges that are stable and preserve their shape over long distances. Nevertheless, following the analysis of Shvidchenko and Pender (2001) and Sukhodolov and Nikora (2012), the present domain is considered to be well suited to cover most of the spectrum of turbulent eddies, which scale with the flow depth H. The relative submergence of the bed particles is H /D = 9. The bulk Reynolds number is the bulk velocity of the flow (Table 1). Except if stated differently, all variables are normalized with H and U b . The friction velocity u τ is defined by extrapolating the linear profile of the total fluid shear stress of the unladen flow down to y = 0, which corresponds to the crest of the immobile layer of particles at the bottom of the domain. This value was used as a reference for all particle-laden flow scenarios in the estimates of the Shields number as explained below. The resulting particle Reynolds number D + = u τ D/ν f is 19.2, indicating that the present simulations fall into the transitionally rough roughness regime. To resolve the viscous length scale at the particle surface throughout the flow, an equidistant Cartesian grid of N x × N y × N z = 2400 × 223 × 1206 with constant isotropic cell size over the whole domain was employed, yielding a total amount of 645 million grid cells. This results in 22.2 grid points per particle diameter and a resolution of + x = u τ x /ν f = 0.86 in terms of reference wall units. Figure 2b illustrates that this resolution, indeed, is fine enough by reporting an instantaneous snapshot of the streamwise velocity, which is very smooth on this grid. A total number N l = 1552 of Lagrangian marker points on the surface of a single particle was used. The fairly large spatial extent of the domain and the high resolution employed required significant computational time. Due to this limitation the Reynolds number had to be kept small, at a value barely above the turbulence threshold for flows over a smooth wall (e.g. Kim, Moin, & Moser, 1987). It is well known that at this Reynolds number both viscous and turbulent effects play an important role in the fluid-particle interaction (Nezu & Nakagawa, 1993). The simulations were performed using 256 cores of a SGI Altix facility and consumed in total more than 600,000 CPU hours. The most important numerical parameters of the simulations are provided in Table 2.

Simulation scenarios and key parameters
To simulate a flow within and above a mobile granular bed, 2000 particles of the same diameter D as those of the fixed sediment bed were released in the flow. Initially, they settle on the bottom and then are transported by the flow as bed-load. Two different scenarios were studied: (i) particles with a larger relative submerged density ρ = (ρ p − ρ f )/ρ f (case HP, "heavy particles"); and (ii) particles with a lower density (case LP, "light particles"), with the defining parameters shown in Table 3. Thus, according to the values of ρ , the case LP exhibits a higher bed mobility, while the case HP shows a lower bed mobility (with both scenarios corresponding to the mobile-bed conditions). Since the amount of mobile particles is equal to 30% of a single layer of the fixed bed, particle protrusion can locally become relatively large, substantially enhancing the actual particle mobility (Fenton & Abott, 1977). The Shields parameter Sh = u 2 τ / (ρ gD) does not account for high protrusion effects and thus serves only as an indicative parameter, rather than an exact criterion of incipient motion. Its estimate for heavy particles (HP) was 0.024 while for the light particles (LP) the Shields parameter was 0.075 determined with the clear-water value of the shear velocity.
The simulations were first run until the erosion and deposition rates were in equilibrium, which was verified using the method presented by Vowinckel, Kempe, and Fröhlich (2013). The corresponding time period is denoted t init . After the equilibrium was achieved, the data for the statistical analysis were stored in terms of entire flow fields over a duration T sample . The corresponding values are provided in Table 3 as multiples of T b = H /U b . A preliminary analysis of these two scenarios in terms of one-point statistics has been presented in  and . In the remainder of this section, a description of the key flow features is given, before the DAM approach is applied to the data in the subsequent sections.

Scenario with heavy particlesa
The heavy particles in scenario HP have a strong tendency to form streamwise clusters of resting particles as visible in Fig. 3. Three characteristic large-scale clusters of particles with an average spacing of 2 H can be observed in the snapshots, albeit with different coherence. While the structure at z ≈ 2 H occupies the entire streamwise extent of the channel, a limitedsize cluster of resting particles with moving particles at its upstream part emerges at z ≈ 0 H. At z ≈ 4 H, the pattern is Table 3 Overview over parameters distinct for the two scenarios addressed in the present study rather irregular, with somewhat scattered particle clusters. The observed particle clusters clearly introduce spatial heterogeneity of the sediment bed in both the streamwise and spanwise directions. The cluster propagation speed differs by orders of magnitude from the fluid velocities, yielding a significant separation between the morphological time scale and the turbulent time scale. Figure 3 illustrates that in response to the particle clusters spatial scales of coherent fluid structures span from the particle scale to the full streamwise and wall-normal extent.

Scenario with light particles
Due to their much higher mobility, the light particles in the scenario LP travel over the fixed bed with a significant streamwise velocity, as displayed in Fig. 4 by their colouring. During the motion, the particles undergo frequent collisions with the bed and other mobile particles. It has been shown in  that the mobile particles form small-scale clusters which are not stable in time. Such a behaviour results in a fairly even distribution of the particles and no clusters of resting particles are observed. Despite the randomness of the distribution of the particles, the increase of large-scale coherent fluid structures is observed compared to the situation of a clear-water flow without mobile particles. This increase, however, is not as noticeable as in the scenario HP.

Background
As mentioned in the introduction, large amounts of data produced by DNS and also by LES require a sound framework that allows data "reduction" to produce a manageable number of interpretable quantities to characterize the flow. For smooth-bed flows researchers and engineers may successfully apply the Reynolds-averaged Navier-Stokes (RANS) framework, which deals with time-or ensemble-averaged variables, but involves no spatial averaging. Although these equations may still be used, in principle, for fixed rough-bed flows, their practical application is problematic due to complex boundary conditions. This inconvenience is removed by spatial averaging of RANS equations that produces double-averaged (in time and space) hydrodynamic equations (e.g. Finnigan, 2000;Giménez-Curto & Lera, 1996;Nikora, Goring, McEwan, & Griffiths, 2001;Raupach & Shaw, 1982). This framework is known as the DAM.
For a more general case of mobile rough-bed flows the conventional Reynolds-averaged hydrodynamic equations are not applicable and thus need to be replaced with a more general form of the double-averaged hydrodynamic equations accounting for the effects of mobile boundaries. A detailed discussion of these issues is provided in Nikora, McEwan, et al. (2007), Nikora, McLean, et al. (2007) and Nikora et al. (2013).

Application of the DAM to the Navier-Stokes equations
To interpret the present DNS data, the double-averaged momentum equation for mobile bed conditions is employed (Nikora et al., 2013). Full details of the derivation as well as considerations of alternative averaging procedures are given in this reference. Here, a short summary of the concept is provided in order to keep the paper self-contained. The consecutive time-space averaging operator for a fluid quantity θ can be defined as: where an overbar denotes time averaging and angular brackets indicate spatial averaging, γ is the clipping (or indicator) function equal to 1 if a point is occupied by fluid and equal to 0 otherwise (e.g. Gray & Lee, 1977), V 0 is the total volume of the averaging domain, T 0 is the total averaging time, φ T = T f /T 0 is the time porosity, and T f is the total time within T 0 during which a given location is occupied by fluid. The parameter φ Vm = V m /V 0 is the space porosity where V m is the volume within V 0 occupied by fluid (even briefly) within the total averaging time T 0 . The symbol with overbarθ and no sub-and superscripts denotes the so-called intrinsic time averaging over temporal sub-domains during which a point in space was actually occupied by fluid, while θ denotes the corresponding intrinsic average in space. The superscript s on the left-hand-side of Eq. (6) denotes the superficial average over the entire volume V 0 and/or the entire averaging time T 0 . As a consequence, the total porosity of the double-averaged mobile granular bed can be defined as φ VT = φ Vm φ T (Nikora et al., 2013).
The averaging operator defined by Eq. (6) gives rise to a modified Reynolds decomposition for intrinsic quantities, since deviations in both time and space are possible, i.e.: In Eq. (7), the prime indicates a deviation of the instantaneous value from its time-averaged value and the tilde indicates deviation of the time-averaged value from its spatially averaged value. To double-average the Navier-Stokes equations, Eqs (6) and (7) need to be supplemented with the so-called double-averaging theorems that link the double-averaged temporal and spatial derivatives with the corresponding derivatives of the double-averaged quantities (Nikora et al., 2013). These theorems involve surface integrals taken over the interface S int between the fluid and mobile granular bed. The combined application of Eqs (6) and (7), and the double-averaging theorems for temporal and spatial derivatives lead to the following double-averaged momentum equation for mobile-bed flows: which constitutes the basis of the analyses presented below. Terms 1 and 2 in Eq. (8) represent local and convective accelerations, respectively. The third term is the momentum supply by the driving volume force, the fourth term results from the pressure gradient, and term 7 accounts for the double-averaged viscous stresses. The fifth and sixth terms are contributions from turbulent and form-induced stresses that originate from the nonlinear convection term in Eq.
(1). The terms described so far are similar to those of the RANS equations (e.g. Rodi, 1993), or for the phase-averaged Navier-Stokes equations (Reynolds & Hussain, 1972). In addition, the eighth and ninth terms represent momentum fluxes, i.e. stresses, due to potential spatial correlations between the local time porosity and time-averaged velocities and emerge from the nonlinear convection term as well. The final two terms in Eq. (8), i.e. the 10th and 11th, are the interfacial terms for pressure and viscous drag on the phase boundary S int arising from the double-averaging theorems, where n i is the unit vector normal to the particle surface and directed into the fluid. Observe that Eq. (8) is an equation for the fluid phase alone. Collision forces act between the particles directly when they touch and hence do not appear in Eq. (8). If the motion of particles is affected by collisions, though, this affects the surface terms and enters in Eq. (8) via this path. The double-averaged motion equation similar to Eq. (8) can also be formulated for the solid phase. Both equations, for fluid and solid phase, are then coupled through interface terms representing physical interactions between phases. An example of such an equation for the sediment conservation following the same philosophy was proposed by Coleman and Nikora (2009). Summing the equations for fluid and sediment motions leads to the single equation describing the motion of the fluid-sediment mixture where the interfacial terms cancel each other and thus vanish.
So far, the dimensions V 0 and T 0 of the double-averaging domain in space and time have not been specified. They cannot be selected completely arbitrarily, without affecting the algebra leading to Eq. (8). Nevertheless, there still is substantial freedom in how to choose these quantities at best for a given flow configuration, which will be discussed in the subsequent Sections 5.2 and 5.3.

Numerical issues
The data from the DNS are available in the form of threedimensional fields of velocity and pressure values at all points of the simulation grid together with all particle centres at the respective instants in time. The computational grid employed by the simulation code is staggered as described in Section 2.1 and illustrated in Fig. 5. The step size of the grid is constant and is the same in all directions. A staggered grid, however, is highly inconvenient for post-processing. Therefore all data are linearly interpolated to the cell-centred points of the grid. The next step is to determine the clipping function γ , also in the form of cell-centre values. The strategy employed is to set γ to zero whenever one of the six points on the faces of this cell is located inside of a solid particle. This approach increases the solid domain by a tiny amount, but was used intentionally to eliminate potential numerical issues which might arise from the IBM technique in the close vicinity of the interfaces. The related increase of the solid domain is so small that its effect is negligible in the averaging process applied later on.
The averaging extent in streamwise direction L 0,x was chosen to be an integer multiple of the step size of the DNS grid x , and an integer fraction of the domain size, L x , such that K L 0,x = L x for some integer K. To obtain a coarser representation, the averaged quantities are evaluated only at certain points X i = L 0,x /2 + (i − 1) L 0,x , i = 1, . . . ,K. All spatially averaged data are then stored on this coarser grid. Due to the Figure 5 Identification of the clipping function γ (x ij , t) in the vicinity of a phase boundary (dashed line) in two dimensions. Arrows represent velocity vector components, red dots correspond to γ (x ij , t) = 0, green dots show γ (x ij , t) = 1, where x ij identifies a cell centre choice of L 0,x , this grid is regular and the averaged data fulfil the same periodic boundary conditions as the original DNS data. Derivatives are then evaluated by applying central finitedifference formulas of second order on this coarser grid. The same procedure was employed to determine averages in y, z and t.

Selection of averaging time
The scenarios HP and LP had run until stationary and uniform flow conditions were reached to minimize the impact of the initial conditions on the statistical results, before data recording for the actual statistical analysis was started in a subsequent second phase. The duration of this second phase needs to be long enough to provide appropriate statistical convergence when computing double-averaged quantities such as the terms in Eq. (8). Since the exact results cannot be obtained a priori, first-and second-order velocity statistics,ū, u u , and u v , were computed for different averaging times T 0 to determine an appropriate duration of the simulation. For sufficiently large T 0 , the statistical parameters approach their expected values and thus the results become fairly independent of T 0 . This condition was tested for three different wall-normal coordinates and six spanwise coordinates at x = 0, which gives a total of 18 test locations distributed evenly over the cross section. This analysis is illustrated in Fig. 6 for two selected spanwise coordinates of scenario HP: (i) z = 2H, corresponding to the ridge extending over the full streamwise length of the simulation domain and stable during the whole simulation time period; and (ii) z = 6H, corresponding to the particle cluster that moved slowly over the fixed bed in the streamwise direction as confirmed by animations of the simulation data. While the statistics computed for the locations above the stable ridge at z = 2H are well converged for all coordinates and statistical quantities explored, the situation is different for the locations above the unstable ridge at z = 6H. At points P 5 and P 6 in the outer flow region, well above moving particles, the first-order velocity statistics approach constant values (Fig. 6a). In the near-wall region, however, the time-averaged velocityū at P 4 starts to increase at T 0 > 100 T b . The situation is even more pronounced for the second-order statistics in Fig. 6b (dotted line of P 4 ). The reason for this sudden transition from low to high values in the statistical data is the propagation of the particle clusters. The visually inspected time evolution of the particle clusters (not shown here) reveal that the position P 4 = (0, 0.11H, 6H ) is significantly influenced by the slow-moving particle cluster until T 0 = 95 T b . This cluster substantially slows down the fluid at P 4 due to the exchange of slow-moving fluid within the cluster and the flow region above. As soon as the cluster has passed this location, the hydrodynamic mechanisms slowing down the fluid in the vicinity of the particles break down and the local fluid velocity and the turbulent fluctuations become stronger. This feature is most pronounced for position P 4 , while for the other points investigated no such drastic changes are visible. These tests illustrate that different regions of the simulation domain can exhibit different statistical behaviour, influenced by either particle clusters or by "clear-water" flow, if particles are only partially mobilized. For the scenario with light particles (LP), the homogeneous distribution of the particles leads to much faster statistical convergence compared to the HP case. Hence, a shorter duration of averaging could be used for the LP case (Table 3).
Following arguments of Nikora et al. (2013), the averaging time must well exceed the turbulent time scale T b , but still has to be much shorter than the morphological time scale of the moving particle clusters. Based on the investigations reported, it was decided to use the total simulation duration as the averaging time for both scenarios to assure that statistical convergence is maximized while the local acceleration term in Eq. (8) is minimized. Hence, the present data are well suited to address a wide range of flow conditions, such as flows over stable roughness features (case HP, z/H = 2), unstable roughness features (case HP, z/H = 4), and homogeneous particle distribution (case LP).

Selection of averaging domain
The spatial averaging domain of volume V 0 should be selected accounting for two requirements. On the one hand, it needs to be fine enough to resolve large-scale spatial heterogeneity (e.g. due to particle clusters). On the other hand, it must be large enough to smooth heterogeneities at the grain scale. This section investigates this issue by considering a range of options for the averaging domain with various sizes.
Due to the boundary conditions employed, an open-channel flow is anisotropic and inhomogeneous in wall-normal direction and thus profiles of turbulence quantities along this direction exhibit strong gradients, particularly near the bed. To properly approximate these gradients, the wall-normal extent of the averaging domain was chosen to be equal to the cell size of the numerical grid in the vertical direction, i.e. L 0,y = y = D/22.2 = 0.005 H, which is the finest resolution possible, even satisfying the prerequisite of a DNS to resolve the smallest scales of the flow and actually involving no averaging in post-processing.
For the horizontal directions, a distribution of the time porosity along the flow is close to homogeneous if the averaging time is sufficient (Fig. 7) while in the spanwise direction the data reveal heterogeneity due to the particle clusters described in Section 3. This spanwise heterogeneity is more pronounced for scenario HP as the particle clusters are formed by mostly resting particles. In scenario LP, the spanwise distribution of the time porosity is much less heterogeneous. The spanwise extent of the averaging domain, hence, needs to account for potential spanwise heterogeneity. Since the spanwise spacing of the particle clusters of scenario HP is around 2 H, an extent of the averaging domain in the spanwise direction larger than 2H = 18 D would smooth the flow heterogeneity due to the particle clusters.
The effect of the size of the averaging domain on the averaged quantities, described in Section 4.2, was explored by calculating wall-normal profiles of flow quantities for different choices of the horizontal extent of the averaging domain. This was done for the spanwise positions z/H = 2, 4, 6 using an averaging domain with a quadratic horizontal shape V q,m , where the subscript m is an integer number identifying the size in multiples of D and q means quadratic, as detailed in Table 4 Journal of Hydraulic Research Vol. 55, No. 2 (2017) Momentum balance in flows over mobile granular beds 199 12H × y × 2D ("Quadratic shape"). Figure 8 shows the wall-normal profiles for different choices of V q,m at z/H = 6, because it was found in Section 5.2 to be the most "delicate" location. In order to describe quasi-steady conditions, the sensitivity to the averaging time at this location can be counterbalanced by a suitable size of the averaging domain. The six options shown in Fig. 8 can be grouped into three pairs: (i) with an averaging domain comparable in size to the particle diameter D, V q,1 and V q,2 ; (ii) with an intermediate extent larger than D but smaller than the flow depth H, V q,4 and V q,6 ; and (iii) with a large size equal to or exceeding the flow depth H, V q,9 and V q,18 . The smallest domains yield the lowest porosity and the strongest fluctuations in the near-wall region, while the large domains exhibit the highest porosity and the lowest fluctuations ( Fig. 8a and  8c). This difference reflects the fact that the larger averaging domains include both ridges and troughs of particle clusters. Hence, an averaging domain with an extent of less than four particle diameters is necessary to properly resolve the depth-scale heterogeneity in the spanwise direction. Option V q,1 , however, shows some waviness in the profile of the velocity variance in the outer flow (Fig. 8c), indicating that this averaging domain is too small to provide adequate statistics for the present analysis. Based on the above analysis, a spanwise averaging over 2 D was selected to smooth the grain-scale heterogeneity while preserving ridge-scale variability. Below we show that the convergence for this spanwise extent can be further improved by increasing the streamwise extent of the averaging domain.
The spanwise extent being fixed, the streamwise extent of the averaging domain was explored using an averaging domain with an elongated horizontal shape V x ,m for the options listed in Table  4, where the subscript m is a number identifying the size in multiples of D and x means averaging along the flow. As with the spanwise averaging, the extent of the averaging domain in the 200 B. Vowinckel et al. Journal of Hydraulic Research Vol. 55, No. 2 (2017)  x direction depends on the specific quantities under consideration. Due to the periodic boundary conditions in the streamwise direction, the averaging domains up to the extent of the computational domain have been explored. A specific pattern to be studied in the present case is the one of streamwise ridges extending over the entire stretch L x , or at least a substantial portion of it. Therefore, it seems that streamwise averaging over the entire extent of the simulation domain is appropriate. To systematically investigate the impact of this issue which is of practical interest for experimentalists, averaging in the spanwise direction over 2 D was retained from the above discussion and averaging in the x direction varied as indicated in Table 4 ("Elongated shape").
Since the location of x = 0 and z = 6 H was identified to be the most sensitive in the computational domain in terms of statistical convergence, the wall-normal profiles at this location are shown in Fig. 9. The statistics of the total porosity and the averaged streamwise fluid velocity converge well for an averaging domain longer than 72 D = 6H = 0.5 L x . The velocity variance is suitably converged in the outer flow (Fig. 9c), but slight deviations of ≤ 15% remain in the near-wall regions, when comparing the results of V x ,72 with the reference domain V x ,108 . At other locations, the differences were substantially smaller. Averaging in x over L x is therefore retained for further analyses as it leads to the largest number of samples.
According to Nikora et al. (2013), a suitable averaging domain must be chosen such that its extent is larger than the "homogenized" bed elements but smaller than larger-scale features one desires to resolve. In the present study, two spanwise scales exist that can be homogenized: (1) the particle scale (with ridges resolved); and (2) the ridge scale (with no largerscale features to resolve). Hence, two averaging domains are employed here: a "local" domain (larger than D but smaller than the spanwise ridge spacing) and a "global" domain (larger than ridge spacing). First, a local averaging domain of size: is chosen to investigate ridge-induced flow heterogeneity in the spanwise direction. In the following, this averaging will be referred to as "local averaging". Second, the total spanwise V 0,1 = L x y L z (10) was used to explore the "globally-averaged" flow properties. This averaging will be referred to as "global averaging". Note that with both of these choices the derivatives in the streamwise direction vanish due to the periodic boundary conditions in the x-direction and the averaging over L x .
6 Momentum balance within the double-averaging framework

Significance of the terms of the double-averaged momentum balance
The aim of the present subsection is to quantify all the terms of the double-averaged momentum equation (Eq. (8)) and then to identify the most significant contributions and those that can be neglected, allowing a simplification of Eq. (8). This should help in deciding on an appropriate form of the momentum equation and the terms of the double-averaged momentum balance that have to be modelled. The required information can be obtained from the present DNS data by computing each individual term in Eq. (8), employing the averaging domain defined by Eq. (9) and evaluating the relative magnitudes of the obtained values. However, this is not possible for terms 10 and 11, which contain integrals over the phase boundary. The IBM yields a smoothed velocity field and does not evaluate the integrals for the viscous force and the pressure force separately (Uhlmann, 2005). Therefore, the data on these local forces are not available and terms 10 and 11 cannot be computed a posteriori for dense particle packings. Term 1 in Eq. (8) drops out as the considered double-averaged flow is steady and therefore the time derivative in this equation vanishes. Terms 2-9 in Eq. (8) were computed employing the spatial averaging domains of Eqs (9) and (10) and the discretization described in Section 5.3. Since the averaging is done over the full length of the computational domain L x , the terms containing streamwise derivatives (i.e. along x 1 ) equate to zero and are thus omitted. The same applies to all terms with j = 1. Tables 5 and 6 assemble the maxima of the absolute values of all computed non-vanishing terms for scenarios HP and LP to provide information on their magnitude and their significance for the momentum balance. To maintain the non-dimensional presentation, the normalization by U 2 b /H is used in the following considerations but not shown explicitly for brevity.
In both scenarios, the momentum exchange between the flow regions within the sedimentary bed and above it is governed by the off-diagonal terms with i = 1 and j = 2. In addition to this, significant convective and form-induced stresses are noted for i = 1 and j = 3 in Tables 5 and 6, indicating a transverse transfer of momentum in the flow.

Spatial variability of the terms of the longitudinal double-averaged momentum balance
The features highlighted in Section 6.1 above can now be discussed in more detail using the locally averaged quantities of the longitudinal momentum balance (i = 1) summarized in Figs 10 and 11. These figures are organized as follows. Quantities related to scenario HP are placed on the left-hand side, while the results for scenario LP are on the right-hand side to allow easy comparison. The spatial distribution of each computed term of Eq. (8) is displayed as a two-dimensional plot of their locally averaged values obtained with an averaging domain defined by Eq. (9). The colour scale is the same for all plots assembled within one figure. A one-dimensional plot on the left of each figure represents the global average obtained with an averaging domain defined by Eq. (10). Figure 10 shows the double-averaged convective acceleration, which is a measure of the transport of momentum by advection. For scenario HP, the convective acceleration (term 2 with i = 1, j = 2) is strongest in the near wall region and at the upper boundary (Fig. 10a). This is most pronounced above the location of a stable ridge at z/H = 2. The one-dimensional profile of the globally averaged quantities shows high absolute values in the bed-load layer 0 < y < 0.2 H. The contour plot shows a pattern of alternating positive and negative values of the convective acceleration in the near-wall region as well as in the outer flow indicating the existence of cells of secondary currents induced by the ridge-type particle clusters on the bed. In the regions of clear-water flow, where the influence of particle-fluid interactions becomes negligible, the advection term diminishes. A similar picture is observed for scenario LP (Fig. 10b), although the convective acceleration does not show such a distinct pattern as in the HP case. The locally-averaged acceleration term with i = 1, j = 3 for both HP and LP does not reveal specific patterns as noted for the case with i = 1 and j = 2. Its globally-averaged values are negligible in both cases HP and LP, as one would expect.
Due to its linear dependency on bed porosity, the momentum supply (term 3, Fig. 11a and 11b) directly reflects the distribution of φ VT reported in Vowinckel, Nikora, et al. (2017). In the outer flow field the momentum supply is constant (as the porosity is equal to 1), while it decreases in the near-wall region, reflecting the decrease in porosity (Fig. 11). The heterogeneity of the momentum supply in the near-wall region is strong for scenario HP (Fig. 11a) and fairly weak for scenario LP (Fig. 11b).
The contributions from turbulent shear stress (term 5) in the clear-water regions are fairly constant and homogeneous for both scenarios (Fig. 11c and 11d), as expected. In the near-wall region, the effects of the ridges become obvious for scenario HP. Here, very low (negative) values of the rate of change of turbulent stress occur around z/H = 1 and z/H = 3 at y = 0, which are the areas surrounding the stable ridge. Such a heterogeneous pattern is not present for scenario LP, where the particle mobility is much higher. The vertical profile of term 5 for the globally averaged quantity shows two distinct peaks at y/H = 0 and y/H = 0.1 for the case HP reflecting the geometry of the fixed bed and a first layer of potentially mobile particles either resting or rolling on the fixed bed. The same quantity for the case LP also shows two peaks, but less pronounced, in particular the upper one.
The two-dimensional plots of the local form-induced stress terms depicted in Fig. 11e and 11f show near-zero values in the outer (clear-water) regions and non-zero (but low) values in the near-wall region for both scenarios. The level of heterogeneity of this term is much higher in the case HP. For this scenario, clusters of resting particles are formed, inducing heterogeneity, which is reflected in the heterogeneity of the rate of change of the form-induced stress term. In the transitional region between the mobile bed and the clear-water layer (0.1 > y/H > 0.2), regions of positive values for the form-induced stress term are observed (Fig. 11e), while their magnitude becomes negligible in the outer flow. The profile of the globally averaged values shows even stronger stress term than obtained from the local averaging, Eq. (9). This effect is due to the larger-scale heterogeneity introduced by bed ridges. For scenario LP, fairly significant values of the stress term are visible at the transition from the fixed bed to the fluid at y ≈ 0 only, for both local and global averaging, with a homogeneous distribution of the local values across the flow.
The values of the viscous term for scenario HP are non-zero at two elevations, where term 7 is sizable (Fig. 11g): one at the upper boundary of the fixed bed and one at y = 0.11 H = 1 D, due to the effect of clusters of resting particles. In the case LP, the viscous term is significant only within the narrow layer just above the fixed bed particles. These distribution patterns could be expected, as the no-slip condition generates the largest gradients in the flow field at the particle surfaces so that only in the near-particle regions term 7 can be significant ( Fig. 11g  and 11h). This term is the one with the highest magnitudes near the boundary among the quantities investigated, which is in line with the studies of turbulent channel flows over smooth walls (e.g. Hoyas & Jiménez, 2006;Kawamura, Abe, & Matsuo, 1999;Kim et al., 1987). The heterogeneity of the viscous term in the spanwise direction is evident for the scenario HP while the scenario LP exhibits homogeneous distributions across the flow due to the high particle mobility.
The present analysis illustrates that similar to the clear-water flow, the momentum balance is governed by momentum supply which is mainly balanced by the effects of turbulent fluctuations in the outer flow and the effects of the viscous stresses in the near-wall region, if global averaging is performed. However, when looking at the two-dimensional plots of the locally averaged quantities, it becomes clear that the bed structures (i.e. particle clusters) introduce heterogeneity in the spatial distributions of flow quantities. When bed forms are present, the convective acceleration and form-induced stress terms play an important role on the local scale. The spatial heterogeneity is also visible in the distributions of the momentum supply and rates of change of turbulent and viscous stresses. The reported observations may be critical for hydrodynamic modelling as one has to account for additional stress terms whenever information about the spanwise heterogeneity is needed. If, however, one is interested in bulk flow behaviour, the influence of the additional stress terms introduced in the double-averaged momentum balance were found to be not as significant.

Discussion and conclusions
Inertial, cohesionless particles of finite size forming a granular bed create a variety of multi-scale spatial heterogeneities in the bed morphology that may significantly affect the flow structure and particle-fluid interactions. Hence, interrelations between flow and mobile particles have to be studied within a framework coupling flow and particles in a rigorous way, which is also flexible enough to account for the individually transported particles as well as for particle clusters of various kinds. In this paper, the framework of the DAM was employed to provide a full description of a flow with mobile, granular beds by analysing data from particle-resolving DNS.
An important issue when applying the DAM is the selection of time and spatial averaging domains. In relation to time averaging it should be highlighted that in the scenario with heavy particles two distinct time scales could be identified: morphological and turbulent time scales. These scales appeared to be well separated, with the former being orders of magnitude larger than the later. This allowed time averaging to be performed over all turbulence and particle scales, preserving at the same time long-term evolution of the flow field due to morphological changes. The situation with the scenario LP was easier as the averaging time well exceeded both turbulent and particle time scales. Furthermore, the flow features also manifest different length scales, covering the grain-scale, the large scale turbulent eddies, morphological forms, and the full extent of the channel. This issue has proved to be very important when choosing the size and shape of the spatial averaging domain. It was shown in detail that some of these length scales introduce heterogeneity in spanwise direction, which is reflected in the double-averaged momentum balance on the local scale. The terms of the momentum balance were thus calculated employing two spatial averaging domains: local averaging domain and global averaging domain.
A full description of the momentum balance was provided by computing all individual stress terms related to the fluid. The analysis shows that form-induced stresses can contribute significantly to the momentum balance and can locally become quite strong. This underlines the necessity for employing an averaging domain that is able to capture the spatial heterogeneities in the flow. Using different averaging domains, it was possible to analyse the effect of clusters on the flow created by mobile and fixed particles. The analysis highlights the capability of the DAM to describe particle-laden flows and introduces a variety of statistical quantities which can be applied to other simulation data of this type. The obtained results may serve as a knowledge base for improving existing flow models for mobile-bed conditions.
In the companion paper by Vowinckel, Nikora, et al. (2017) the knowledge gained from the present study is further exploited to compute momentum fluxes, the bulk momentum balance, and the modified resistance of the mobile granular bed to the flow as well the structure of secondary currents due to ridge effects.
f = ( f x , f y , f z ) T = driving volume force (m s −2 ) f IBM = IBM source term (m s −2 ) g = gravity acceleration (m s −2 ) H = flow depth above fixed particle crests (m) I = identity matrix (-) I p = moment of inertia of a particle (kg m 2 ) L x , L y , L z = extent in streamwise, wall-normal, and spanwise directions (m) M p = moment due to particle interactions (kg m 2 s −2 ) m p = particle mass (