Spatially-averaged flows over mobile rough beds: equations for the second-order velocity moments

ABSTRACT The double-averaging methodology is used in this paper for deriving equations for the second-order velocity moments (i.e. turbulent and dispersive stresses) that emerge in the double-averaged momentum equation for incompressible Newtonian flows over mobile boundaries. The starting point in the derivation is the mass and momentum conservation equations for local (at a point) instantaneous variables that are up-scaled by employing temporal and spatial averaging. First, time-averaged conservation equations for mass, momentum, and turbulent stresses for mobile bed conditions are derived. Then, the double-averaged hydrodynamic equations obtained by spatial averaging the time-averaged equations are proposed. The derived second-order equations can serve as a basis for the construction of simplified mathematical and numerical models and for interpretation of experimental and simulation data when bed mobility is present. Potential applications include complex flow situations such as free-surface flows over vegetated or mobile sedimentary beds and flows through tidal and wind turbine arrays.


Introduction
Overland, river, coastal and atmospheric flows frequently exhibit high levels of multi-scale spatial heterogeneity due to geometrically complex and, occasionally, mobile boundaries. This heterogeneity affects the instantaneous and time-averaged flow fields, especially in the region near the bounding surface. Various averaging techniques such as ensemble, temporal and/or spatial averaging have been typically applied to cope with this flow heterogeneity (e.g. Ishii & Hibiki, 2006;Jackson, 2000;Monin & Yaglom, 1971;Whitaker, 1999). Applied to the original equations of fluid mass and momentum conservation, the averaging enables up-scaling so the newly obtained equations describe somewhat homogenized fields of averaged hydrodynamic variables.
Combined time and spatial averaging (known as doubleaveraging) has been extensively used in the studies of hydraulically rough-bed flows as it explicitly (i) incorporates the effects of form drag and viscous friction in the double-averaged equations, and (ii) accounts for roughness-induced small-scale heterogeneity in turbulent and mean flow fields (e.g. Giménez-Curto & Lera, 1996;Nikora et al., 2007;Nikora, Ballio, Coleman, & Pokrajac, 2013;Nikora, Goring, McEwan, & Griffiths, 2001;Pedras & de Lemos, 2000;Raupach & Shaw, 1982). This type of heterogeneity in the instantaneous and time-averaged flow fields is presented in the double-averaged equations as turbulent stresses and form-induced or dispersive stresses (which are second-order velocity moments, e.g. Giménez-Curto & Lera, 1996;Monin & Yaglom, 1971;Raupach & Shaw, 1982). Equations for these fluid stresses, known as second-order equations, represent an important part of the theoretical background underpinning modern studies of turbulent wall-bounded flows.
Time-averaged second-order equations describing the balances of mean energy and turbulent energy (i.e. normal stresses) were originally introduced in the pioneering work of Reynolds (1895). Later they were expanded to also cover shear stresses and have served as a foundation for data interpretation and a variety of turbulence modelling approaches, such as families of k-models and algebraic stress turbulence models (e.g. Monin & Yaglom, 1971;Nezu & Nakagawa, 1993;Rodi, 1984). In the time-averaged second-order equations, the drag forces are accounted for as the boundary conditions rather than included in the equations themselves. Therefore, such equations become unsuitable, in a practical sense, for understanding flows over complex boundaries such as gravel or dune beds in rivers. This limitation stimulated the development of the doubleaveraged second-order equations where drag forces and their effects appear in the equations explicitly, as a result of spatial averaging (e.g. Brunet, Finnigan, & Raupach, 1994;Katul & Albertson, 1998;Mignot, Barthélemy, & Hurther, 2008, 2009Pedras & de Lemos, 2001;Raupach & Shaw, 1982;Wilson & Shaw, 1977).
Another important feature that emerges in the doubleaveraged equations is form-induced stresses, which are covariances of the deviations of the time-averaged velocities from their double-averaged counterparts, similar to the turbulent stresses which are co-variances of the deviations of the instantaneous velocities from their time-averaged counterparts. Although a few studies highlighting the importance of the form-induced stresses in rough-bed flows have been reported (e.g. Aberle, Koll, & Dittrich, 2008;Coleman, Nikora, McLean, & Schlicke, 2007;Cooper, Aberle, Koll, & Tait, 2013;Manes, Pokrajac, Coceal, & McEwan, 2008;Vowinckel, Nikora, Kempe, & Fröhlich, 2017a, 2017b, little attention has been paid to the balance of the form-induced stresses (normal and shear). Although a form of this balance equation was proposed a few decades ago (Raupach & Shaw, 1982), it has only recently been used for studying roughness effects in open-channel flows (Yuan & Piomelli, 2014).
In general, the double-averaged second-order equations include three types of interlinked equations: (1) balance equations for mean momentum fluxes (or stresses), including mean energy balance, (2) balance equations for spatially-averaged turbulent stresses, including turbulent kinetic energy balance, and (3) balance equations for form-induced stresses, including form-induced kinetic energy balance. The currently available double-averaged second-order equations have been derived for fixed boundaries and, strictly speaking, cannot be applied if flow boundaries are both complex and mobile; thus, their application to questions related to river beds during floods or vegetation canopies that constantly fluctuate as a result of interaction with turbulent eddies, is limited. Indeed, boundary motion can cause the discontinuity in the near-boundary instantaneous fluid properties within the averaging period, and thus conventional time-averaging and corresponding time-averaged and double-averaged equations need to be modified to account for boundary mobility. All three types of equations (for mean, turbulent, and form-induced stresses) are important as they underpin any consideration of energy fluxes in turbulent flows over complex boundaries, and are therefore required to provide a proper theoretical foundation for dealing with mobile-boundary flows.
In response to this need, the current paper presents the derivation of equations describing the balances of the mean, form-induced, and spatially-averaged turbulent stresses and energies for mobile-boundary flows. As a necessary initial step, the time-averaging operation is refined to account for the rough-bed mobility effect and obtain time-averaged hydrodynamic equations for mobile-boundary flows. Spatial averaging is then applied on the time-averaged equations to obtain the second-order equations for stresses that complement the doubleaveraged equations of conservation of fluid mass and momentum proposed by Nikora et al. (2013) for mobile-bed conditions. The general approach of Keller and Friedmann (1924; see also Monin & Yaglom, 1971) is followed for the derivation of the equations for the second-order velocity moments for both time-averaged and double-averaged equations.
The paper is structured as follows. Definitions of the averaging operations and the associated quantities of bed porosity (or roughness geometry functions) used in the derivations are presented in Section 2. The local (at a point) time-averaged hydrodynamic equations are given in Section 3, as these are required to obtain the time-and spatially-averaged hydrodynamic equations that are presented in Section 4. An example to illustrate the potential application of the proposed equations, under appropriate simplifying assumptions, is given in Section 5. Conclusions are presented in Section 6. Appendix 1 provides supplementary details on the conditions and theorems that are necessary for averaging operations employed in the derivation process outlined in Appendix 2.

Background
This Section briefly reviews key definitions of the time-and double-averaging operations for mobile-boundary flows that are used for the derivation of the equations presented in Sections 3 and 4. Additional information on the averaging conditions and theorems is given in Appendix 1, while details regarding the double-averaging methodology as followed here can be found in Nikora et al. (2007Nikora et al. ( , 2013.

Averaging operators: definitions
We first define the operation of time averaging at a point x i . The superficial and intrinsic time averages of a random function θ (x i , t) are defined, respectively, as (Nikora et al., 2013): where θ is velocity, pressure or any other hydrodynamic variable, overbar denotes the time-averaging operation, superscript s denotes superficial averaging, x i and t are space and time coordinates, τ is a "local" time coordinate used in the integrand, γ is an indicator of the point occupancy by the fluid (i.e. γ = 1 when the domain is occupied by the fluid and γ = 0 otherwise), T f refers to the sum of periods during which a point x i is occupied by fluid, and T 0 is the total averaging time (T f ≤ T 0 ). The superficial and intrinsic time averages are linked through the local time porosity φ T = T f / T 0 (Nikora et al., 2013), as expressed by Eq. (1c). Equations (1a)-(1c) enable the operation of time averaging over a time period T 0 even if this period includes time intervals or instants when the point x i had not been occupied by the fluid. For brevity, the variable dependence on the space and/or time coordinates is not shown hereafter. Applying the operators of superficial and intrinsic spatial averaging (Gray & Lee, 1977) on the first and second operators in Eqs (1a) and (1b), respectively, yields: where the angular brackets denote spatial averaging, V 0 is the total volume of the spatial averaging domain that in general can be subdivided in two subdomains: (1) a subdomain of volume V m every point of which has been visited by the fluid, at least once, within the period T 0 , and (2) a subdomain of volume (V 0 − V m ) that has not been visited by the fluid within the total averaging time T 0 . The expressions in Eq.
(2) correspond to the definitions of superficial and intrinsic double averages given in Nikora et al. (2013). The superficial and intrinsic double averages are linked through a relation θ s s = φ Vm φ Tθ , where φ Vm = V m / V 0 denotes the spatial porosity based on the nonzeroγ s . The space-time porosity is defined as φ VT = φ Vm φ T (Nikora et al., 2013).
When applied to the hydrodynamic equations, the operators of Eqs (1) and (2) need to be supplemented by the averaging conditions (known as the Reynolds conditions or rules) and theorems detailed in Appendix 1. The application of the Reynolds conditions requires the selection of an appropriate averaging time and spatial domain. The averaging time T 0 should be much larger than the characteristic time scale of turbulent fluctuations θ = θ −θ but much smaller than the characteristic time scale ofθ. For example, in rivers T 0 would be much larger than the ratio of the flow depth to the flow velocity (i.e. prevailing turbulent scale) but much smaller than the duration of a flood. Thus, a sufficiently wide separation between these two time scales is required to secure the applicability of Reynolds conditions. An appropriate shape of the spatial averaging domain is a thin slab parallel to the mean bed, with a vertical size much smaller than the roughness height. The bed-parallel dimensions of the averaging domain should appreciably exceed the roughness length scales along and across the flow but be sufficiently smaller than channel-scale of bed level fluctuations. Applying, e.g., this rule for a gravel-bed river we would define an averaging domain as being much thinner than a bed particle diameter d, with bed-parallel dimensions being a few times larger than d but much smaller than a channel width. Similar to time averaging, a wide separation between characteristic roughness scales and the scale of variation of the spatially-averaged value is required for spatial averaging (Eq. (2)) to apply.

Decomposition of instantaneous and time-averaged variables
Another element for the derivation of the time-and doubleaveraged second-order equations, in addition to Eqs (1) and (2), is the decomposition of the instantaneous variables θ into meanθ and fluctuating θ quantities, known as the Reynolds decomposition, i.e. θ = θ −θ. Application of the Reynolds decomposition to the time-averaged momentum flux φ T u i u k ((i) in Eq. (3a)) results in the following expansion: where quantities (ii) and (iii) represent the contributions of the mean and fluctuating velocity fields to the time-averaged momentum flux, respectively. A spatially-averaged analogue of Eq. (3a) is given as: where the meaning of quantities (II) and (III) is similar to terms (ii) and (iii) in Eq. (3a). The relation of Eq. (3a) involves time averaging and no spatial averaging (at least in the framework of the continuum mechanics), while the relation of Eq. (3b) may be viewed as a spatially "upscaled" version of Eq. (3a).
While temporal averaging is sufficient for studying hydraulically smooth-bed flows, consideration of hydraulically roughbed flows requires spatial averaging and decomposition of a time-averaged quantityθ into its spatially-averaged value θ and a spatial fluctuationθ , i.e.θ = θ +θ . The latter relation is known as a modified Reynolds decomposition. Spatial averaging is necessary to cope with the roughness-induced heterogeneity of the time-averaged velocity field, i.e. to make timeaveraged heterogeneous fields "smoother". Using the modified Reynolds decompositionū i = ū i +ũ i , the double-averaged momentum flux φ Vm φ Tūiūk ((II) in Eq. (3b)) may be decomposed as: where quantities (IIA) and (IIB) represent contributions to the total spatially-averaged momentum flux φ Vm φ Tūiūk from the double-mean ( ū i ) and form-induced (ũ i ) velocity fields, while terms in (IIC) reflect potential contribution of the interrelations between the time porosity φ T and the time-averaged velocity. For convenience, the quantities (IIA), (IIB) and (IIC) are referred to as the double-mean, form-induced, and porositycorrelation stresses (which may be either zero or non-zero), respectively. From Eqs (3b) and (4), a decomposition for the double-averaged kinetic energy 1 / 2 φ Vm φ T u i u i follows as: where repeated indices imply summation. The balance equations for the local mean and turbulent stresses (quantities ii and iii in Eq. (3a)) as well as for the double-mean, form-induced and spatially-averaged turbulent stresses (quantities IIA + IIC, IIB and III in Eqs (3b) and (4), respectively) are presented in Sections 3 and 4 below. The key steps of their derivation are outlined in Appendix 2. The general approach of Keller and Friedmann (1924; see also Monin & Yaglom, 1971) is followed, starting from the mass and momentum conservation equations.

First-and second-order time-averaged hydrodynamic equations
The derivation of the time-averaged hydrodynamic equations for a Newtonian fluid starts with consideration of the mass and momentum conservation equations for instantaneous variables (e.g. Whitaker, 1968), i.e.: where ρ is fluid mass density, ν is fluid kinematic viscosity, u i is the ith component of the velocity vector, p is fluid pressure, and g i is the ith component of the gravity acceleration. These conservation equations are averaged over time first to obtain the time-averaged mass and momentum conservation equations (i.e. first-order equations), and then the balance equations for the mean and turbulent momentum fluxes (i.e. second-order equations) are derived. The final equations are reported below, while details of their derivation are provided in Appendices 1 and 2.

Time-averaged mass and momentum conservation equations
Applying the operators of Eq. (1), Reynolds' conditions of Eq. (A1), the averaging theorems (A5), and accounting for the no-slip condition, one can derive the following time-averaged equations: where the Reynolds decomposition u i =ū i + u i is used in the presentation of the time-averaged convective term, δ denotes the one-dimensional analogue of the Dirac delta function, t s is a "transition" time instant when γ at a given point changes from 0 to 1, v i is velocity vector of the fluid-non-fluid interface, and n i is a normal unit vector at the interface (directed into the fluid). More details on the variables involved in Eq. (8) are given in Appendix 1. The first and second terms in Eq. (8) are the mean local and convective accelerations, respectively. The third term represents the body force (gravity), with the fourth, fifth and sixth terms expressing the contributions to the time-averaged momentum balance of turbulent stresses, mean pressure and mean viscous stresses, respectively. The final two terms in Eq. (8) reflect the change of fluid momentum (per unit mass) at the transition time moments (non-fluid to fluid and vice versa) due to the effects of pressure drag and viscous friction at the interface. The time-averaged continuity Eq. (7) can be expressed as: Equation (9) suggests that temporal and spatial variations of time porosity φ T can cause non-zero strain rate ∂ū i /∂x i , differently from the conventional time-averaged continuity equation where this strain rate is identically zero.

Time-averaged conservation equations for mean and turbulent momentum fluxes
The equations for the meanū iūk and turbulent momentum fluxes u i u k can be derived from the instantaneous and time-averaged forms of the mass and momentum conservation equations, as outlined in Appendix 2. The equation for the mean stressū iūk is obtained as: Terms 1 and 2 in Eq. (10) represent the time rate of change of the mean stresses, term 3 is a source term due to gravity, and terms 4, 5 and 6 are transport terms related to turbulent stresses, mean pressure and mean viscous stresses, respectively. Term 7 is conventionally interpreted as turbulence production due to the work of turbulent stresses against mean velocity gradients. Terms 8 and 9 express redistributive (i.e. inter-component transfer of energy) and dissipative effects, respectively. Terms 10 and 11 are source/sink terms associated with the flow interruption by non-fluid inclusions due to boundary mobility. The balance equation for the mean energy 1 / 2ūiūi can be obtained from Eq. (10) by imposing i = k, i.e.: 1 2 ∂φ Tūiūi ∂t Terms 1 and 2 in Eq. (10) represent the local and convective time rate of change of the mean energy; terms 3, 4, 5 and 6 represent the work rates of gravity, turbulent stresses, mean pressure and mean viscous stresses on the mean flow, respectively; terms 7, 8 and 9 represent the work of turbulent stresses, mean pressure and mean viscous stresses against the mean strain-rate, respectively. Conventionally, the seventh term is interpreted as mean energy exchange between mean flow and turbulence while the ninth term describes conversion of mechanical energy into heat. Terms 10 and 11 represent the gain or loss of mean energy at the transition time moments due to the boundary mobility.
The equation for the turbulent stress u i u k can be obtained by subtracting the equation forū iūk from the equation for the total stress u i u k (Appendix 2). The result can be written as: Terms 1 and 2 in Eq. (12) correspond to the time rate of change of turbulent stresses, term 3 is the turbulent stress generation, and terms 4, 5 and 6 represent the transport by the fluctuating velocity, pressure and viscous stresses, respectively. Term 7 describes the stress redistribution by pressure fluctuations, term 8 is responsible for dissipative effects, and terms 9 and 10 are source/sink terms associated with the flow intermittency, i.e. intermittent occupation of a point under consideration by nonfluid matter (e.g. moving gravel particles or fluctuating aquatic plants). The equation for the turbulent energy 1 / 2 u i u i is obtained from Eq. (12) as follows: In Eq. (13), terms 1 and 2 represent the rate of change of turbulent energy, term 3 represents the work of turbulent stresses against the mean strain rate, interpreted as the conversion of the mean energy to turbulence energy (or vice versa), term 4 describes the turbulent convection and terms 5 and 6 represent the mean work rate of the turbulent pressure and viscous stresses on the turbulent fluctuations, respectively. Terms 7 and 8 denote the mean work of turbulent pressure and viscous stresses against the turbulent strain rate. Terms 9 and 10 express the mean gain or loss of the turbulence energy at the transition time moments (i.e. energy transfer to the boundary or vice versa). The intrinsic forms of Eqs (8) and (10)-(13) can be obtained by applying Leibniz's product rule on the left-hand side terms and then dividing the equations by φ T . For brevity, these are not presented in this paper. The effects of moving flow boundaries and the associated flow intermittency are expressed in the time-averaged equations through the time porosity and the two source/sink terms that appear at the end of Eqs (8) and (10)-(13). Note that due to Eq. (9), the work of pressure on mean and turbulent strain-rates, i.e. term 8 in Eq. (11) and term 7 in Eq. (13), are retained, whereas for the case of immobile boundary these two terms are zero. The time-averaged equations for flows over mobile boundaries can be used directly or they can underpin the development of the double-averaged equations, as outlined in Appendix 2 and summarized in the following section.

Double-averaged hydrodynamic equations
Time-averaged hydrodynamic equations involve no spatial averaging and thus relate to local (at a point) quantities. As detailed in the introduction, their application for describing flows over complex and mobile boundaries, like gravel beds during floods, is impractical, so a second averaging step, in the spatial domain, is needed. This step may be viewed as a homogenization procedure that upscales flow variables from a "point scale" where they are highly heterogeneous to a larger scale where their behaviour is much smoother (e.g. Nikora et al., 2013).
The equations for the second-order velocity moments, i.e. quantities (IIA), (IIB) and (III) in Eqs (3b) and (4), are needed in order to identify the key mechanisms affecting the double-mean, form-induced and turbulent stresses and their interrelations with the boundary motion. The derivation of the double-averaged second-order equations utilizes the original mass and momentum conservation equations (Eq. (6)), their time-averaged forms (summarized in Section 3, Eqs (7) and (8)), and double-averaged forms (Nikora et al., 2013). To complement Eqs (6)-(8), the double-averaged mass and momentum conservation equations are summarized in Section 4.1. The second-order doubleaveraged equations (for momentum fluxes) are presented in Section 4.2, with key steps in their derivation outlined in Appendix 2.

Double-averaged continuity and momentum equations
The double-averaged forms of the mass and momentum conservation equations for mobile-bed conditions can be expressed as: where all variables have been defined in the previous sections. Equation (14) is a double-averaged counterpart of the left-side relation in Eq. (6) or, equivalently, a spatially-averaged form of Eq. (7). Terms 1 and 2 in Eq. (15) represent the local and convective parts of the double-averaged acceleration, term 3 represents the body force (gravity), and terms 4, 5 and 6 reflect the contributions of the turbulent, form-induced and mobility-induced stresses, respectively. Terms 7 and 8 express the effects of the double-averaged pressure and viscous stresses, while terms 9 and 10 are the mean pressure drag and viscous friction terms, respectively (Nikora et al., 2013).

Double-mean stresses
The balance equation for the double-mean stress φ VT ū i ū k and the associated porosity correlation contribution φ Vm φ Tũi ū k + φ Vm φ Tũk ū k is derived using the double-averaged mass and momentum conservation equations together with the identity (B5) that expands the partial time derivatives of φ VT ū i ū k and φ Vm φ Tũi ū k + φ Vm φ Tũk ū k (quantities (IIA) and (IIC) in Eq. (5) and Eq. (B5) in Appendix 2), i.e.: Terms 1 and 2 in Eq. (16) represent the time rate of change of the double-mean stresses while terms 3 and 4 quantify the rate of change of the porosity-correlation stresses. Terms 5, 13 and 14 occur in both balance equations, for double-mean stresses (Eq. (16)) and for the form-induced stresses (Eq. (18) below, terms 3, 5, and 6, respectively), but with opposite signs. Thus, these terms express the energy exchanges between these two energy balances. Term 6 represents the effect of gravity. Terms 7, 8 and 9 describe the transport of turbulent, form-induced and mobilityinduced stresses, while terms 10 and 11 are the pressure and viscous transport terms. Term 12 is present in both Eq. (16) and in the balance of the spatially-averaged turbulent stresses (Eq. (20)) with opposite signs. For this reason, term 12 is interpreted as a "bridge" for energy exchanges between mean and turbulent flow fields. Term 15 reflects the effect of the pressure acting on the strain rates, while term 16 quantifies the viscous dissipation. The interfacial source/sinks of mean stresses are expressed by terms 17 and 18. The equation for the double-mean energy 1 / 2 ū i ū i is obtained from Eq. (16) by choosing i = k, i.e.: Terms 1 and 2 in Eq. (17) are the local and convective time rates of change of 1 / 2 ū i ū i , with terms 3 and 4 representing the effects of the local and the convective parts of the porositycorrelation. Term 5 defines an energy exchange with the balance of the form-induced energy (Eq. (20)). Term 6 describes the energy supply from the potential energy of gravity, while terms 7, 8, 10 and 11 represent the work of turbulent and form-induced stresses, double-averaged pressure and double-averaged viscous stresses on the double-averaged velocity field, respectively. Term 9 represents the spatial transfer of the double-mean energy 1 / 2 ū i ū i by the form-induced flow and it is associated with the boundary mobility, i.e. with potentially non-zero porosityvelocity correlation. Terms 12, 13, 14, 15 and 16 relate to the work of the spatially-averaged turbulent, form-induced, and porosity-correlation stresses, double-averaged pressure, and double-averaged viscous stresses on the double-mean strain rate, respectively. Finally, terms 17 and 18 define the work of mean pressure drag and viscous friction at the fluid-non-fluid interfaces on the double-mean flow.

Form-induced stresses
By subtracting Eq. (16) from the equation for the spatiallyaveraged mean stresses φ Vm φ Tūiūk (following Eq. (4) and Eq. (B6) in Appendix 2), the equation for the form-induced stresses can be derived as: Journal of Hydraulic Research Vol. 58, No. 1 (2020) Equations for stresses for mobile-boundary flows 141 Terms 1 and 2 in Eq. (18) represent the rate of change of the form-induced stresses φ Vm φ Tũiũk . Term 3 of Eq. (18) also appears in Eq. (16) as term 5, but with an opposite sign, as mentioned in relation to Eq. (16). For this reason, it is interpreted as a mobility-induced cross-energy exchange between balances in Eqs (18) and (16). Term 3 may not be zero if there is a spatial correlation between a form-induced velocityũ i and time porosity φ T . Likewise, terms 5 and 6 represent the covariance exchange between double-mean and form-induced balances of Eqs (18) and (16). Term 4 shows that the potential energy of gravity may be directly supplied to the form-induced stress balance, as in the double-mean energy balance in Eq. (17) In Eq. (19), terms 1 and 2 describe the local and convective change of the form-induced energy. Term 3 represents a part of form-induced energy exchanged with the double-mean energy balance of Eq. (17), where it is seen as term 5. Term 4 describes the potential supply of energy from gravity, if φ Tũi > 0. Terms 5 and 6 relate to the work of form-induced and porositycorrelation stresses against the bulk strain rate. The convection of the form-induced energy is described by term 7. Terms 8, 9 and 10 reflect the mean work rate of the turbulent stresses, mean pressure, and mean viscous stresses on the form-induced velocity field, respectively. The work of turbulent stresses, mean pressure, and mean viscous stress on the form-induced strain rates are accounted for by terms 11, 12 and 13, respectively. Terms 14 and 15 represent the work of pressure and viscous stresses on the interfacial surface and have the meaning of the source/sink terms in the form-induced energy balance. Note that the sum of Eqs (17) and (19) describe the balance of the spatially-averaged mean kinetic energy 1 / 2 ū iūi .

Spatially-averaged turbulent stresses
Integrating the local stress balance of Eq. (12) over V m , dividing by V 0 , and applying the double-averaging theorems given in Appendix 1 yield the balance for the spatially-averaged turbulent stresses, i.e.: Journal of Hydraulic Research Vol. 58, No. 1 (2020) The interpretation of the terms in Eq. (20) is similar to that of the terms in Eq. (12), although here it relates to the spatiallyaveraged quantities. Terms 1 and 2 represent the rate of change of the spatially-averaged turbulent stresses. Terms 3 and 4 express the cross-energy (covariance) exchange between Eq. (20) and Eqs (16), (18), respectively. Terms 5 and 6 describe the convective transport by the form-induced and turbulent velocity fields. Terms 7 and 8 are pressure and viscous transport terms. Terms 9 and 10 express the redistributive and dissipative effects. The final two terms, 11 and 12, are the source/sink terms associated with the boundary motion. The equation for the spatially-averaged turbulent kinetic energy 1 / 2 φ Vm φ T u i u i is obtained from Eq. (20) by imposing i = k, i.e.: 1 2 The local and convective rates of change of the spatiallyaveraged turbulent energy are given by terms 1 and 2 in Eq. (21). Term 3 represents the work of the spatially-averaged turbulent stresses against the double-mean strain-rate, while term 4 describes the mean work rate of the turbulent stresses against the form-induced strain-rate. Terms 5 and 6 reflect the form-induced and turbulent convection effects. Terms 7 and 8 represent the work rate of turbulent pressure and viscous stresses, while terms 9 and 10 express the work of pressure and viscous stresses against the turbulent strain rate. Terms 11 and 12 are source/sink terms that relate to the work of pressure and viscous stresses on the interfacial (fluid-non-fluid) surface. In other words, these two terms highlight possibilities for the energy exchange between fluid flow and its moving boundary.

Discussion
The second-order equations reported in this paper supplement the double-averaged equations of mass and momentum conservation proposed in Nikora et al. (2013) for the study of mobile-bed flows (e.g. open-channel flows over mobile granular beds or terrestrial flows over wind-turbine arrays). The spatiallyaveraged turbulent stress, form-induced stress, and drag terms which arise in the double-averaged momentum equation due to averaging are the unknown variables that need to be parameterized or modelled. In this paper, we propose the equations describing the dynamics of the spatially-averaged turbulent and form-induced stresses. These equations can serve as the basis for the construction of predictive models and for the interpretation of experimental data in the studies of mobileboundary flows. The interrelations between the double-mean, form-induced and turbulent components, as well as between these and boundary motions, should also be considered and modelled. The roughness mobility effects and the boundary energy fluxes are features that differentiate the proposed equations from the double-averaged equations reported by others. These features are highlighted in the following discussion, along with some examples of applications of the proposed equations (after employing appropriate simplifying assumptions).

Roughness mobility effects
The mobility of the flow boundaries causes an intermittent occupation of the near-bed domain by either the fluid or bed roughness elements. This intermittency is tracked using the clipping function γ (x i , t) that is defined as 1 when a point x i is occupied by the fluid and 0 otherwise. Using function γ we define the "fluid" averaging period T f (x i , t) as the collection of time periods when the fluid occupies a particular point x i of the domain, and then employ the ratio T f / T 0 = φ T as a measure of boundary mobility. As shown by Eq. (9), the variation of φ T in time and space can impose a non-zero mean strain rate (∂ū i /∂x i ) that may be interpreted as a "compressibility" effect in the time-averaged flow field. Due to Eq. (9) and the fluid incompressibility, i.e. ∂u i /∂x i = 0, a non-zero turbulent strain-rate ∂u i /∂x i is implied. Consequently, the associated work of mean and turbulent pressure on the mean and turbulent strain-rates, attributed to terms 8 and 7 in Eqs (11) and (13) respectively, is retained in the balances of mean and turbulent energy. After spatial averaging, similar terms can be noted in the equations of double-mean, form-induced and turbulent energy balances, i.e. term 15 in Eq. (17), term 12 in Eq. (19) and term 9 in Eq. (21). In addition, spatial correlations between time-averaged velocity and porosity φ T arise in the double-averaged equations that may be interpreted as momentum fluxes due to the boundary mobility. Furthermore, the intermittent occupancy of the domain points by fluid within the averaging time period can cause changes in fluid mean momentum (losses or gains), mean energy, and turbulent energy when the domain occupancy changes. These sharp changes are attributed to the final two terms in each of Eqs (10)-(13). The effects of time porosity φ T and the boundary mobilityrelated time integrals constitute the main difference between the proposed time-averaged equations for mobile-boundary flows and the conventional Reynolds-averaged equations for fixedboundary flows. Operators and theorems for time-averaging over variables that may not be continuous within the averaging time period have been proposed earlier by Ishii and Hibiki (2006) and others (e.g. Drew & Passman, 1999). The conceptual difference between their approach and the one proposed in this paper is the use of the distribution function γ that allows the transition from the time-averaged equations for mobileboundary flows to those for fixed-boundary conditions. With the current approach, the problem associated with a possible zero velocity v i of the boundary in the denominator of the interfacial terms in the equations of Drew and Lahey (1979) is avoided, as for γ = 1 its derivatives are ∂γ /∂t = 0 and ∂γ /∂x i = 0 and the gain/loss integral terms do not emerge. For the fixed-bed flows, the time porosity is φ T ≡ 1 (Nikora et al., 2013) and the operations of time-averaging and differentiation commute. Then, the proposed time-averaged equations coincide with the conventional Reynolds-averaged equations. The effects of the roughness geometry and bed mobility are introduced in the doubleaveraged equations for the second-order moments through the roughness geometry functions (i.e. time and spatial porosities) and the work of pressure and viscous drag on the interfacial surfaces. These are the main differences between the proposed double-averaged equations and the equations developed by others for the form-induced energy (Raupach & Shaw, 1982) or for the spatially-averaged turbulent stresses/energy (Brunet et al., 1994;Mignot et al., 2008;Pedras & de Lemos, 2001).

Some physical interpretations
A preliminary interpretation of the terms involved in the doubleaveraged second-order Eqs (17), (19) and (21) can be based on the theoretical definitions of the rate of work and strain power (Malvern, 1969). Given the form of Eqs (17), (19) and (21), kinetic energy is supplied to the flow from the potential energy of gravity and is directed to the double-mean and forminduced energy balances through the right-hand-side terms of the following relation: Due to the work rate of pressure, viscous stresses, and turbulent stresses, mechanical energy is redistributed over space. Through the work of viscous stresses against the strain rate a part of the mechanical energy is converted into heat, removing energy from the double-mean, form-induced and turbulent fields, i.e.: fluid kinetic energy dissipated from the double-mean field (Eq. (17)) Owing to the work of pressure and viscous stresses on the interfacial surface, the energy is exchanged between the three constituting flow fields (double-mean, form-induced, and turbulent) and boundary motions, i.e.: s work of pressure drag and viscous friction on the double-mean flow (Eq. (17)) work of hydrodynamic forces on the form-induced flow (Eq. (19)) The modified Reynolds decomposition gave rise to the terms that do not emerge in the equation for the total double-averaged kinetic energy 1 / 2 φ Vm φ T u i u i (sum of Eqs (17), (19) and (21)). These are the terms 13 and 14 in Eq. (17) and equivalent terms 5 and 6 in Eq. (19), term 12 in Eq. (17) and corresponding term 3 in Eq. (21), and term 11 in Eq. (19) and its counterpart term 4 in Eq. (21). Such terms represent the energy exchange between the respective balances of the double-mean energy, form-induced energy, and turbulent energy. Figure 1 illustrates the energy exchanges between the double-mean, form-induced and turbulent fields, as well as the conversion of the kinetic energy to heat and the energy transferred through the action of mobile boundaries.
The exchange between the double-mean and turbulent flow fields φ Vm φ T u i u j (∂ ū i /∂x j ) (i.e. term 3 in Eq. (21) and term 12 in Eq. (17)) may be interpreted as the production of the turbulent energy, as in considerations of the conventional Reynolds-averaged equations. This interpretation is supported by estimates for fixed-gravel-bed flows (e.g. Mignot, Barthelemy, & Hurther, 2009;Yuan & Piomelli, 2014). The energy exchange between the form-induced and turbulent fields φ Vm φ T u i u j (∂ũ i /∂x j ) (term 4 in Eq. (21) and term 11 in Eq. (19)) are assessed in Yuan and Piomelli (2014) where they are considered as bidirectional energy exchange depending on the characteristic scales present in the flow. In an earlier study, Mignot et al. (2009) noted that this term may be insignificant compared to the turbulence production term 3 in Eq. (21). Atmospheric flow studies (e.g. Raupach & Shaw, Figure 1 Simplified representation of the energy fluxes between the double-mean, form-induced and turbulent flow fields, conversion of the kinetic energy to heat and the energy exchange (source or sink of energy) with the mobile boundaries of the flow (where T ik /ρ = −pδ ik /ρ + ν∂u i /∂x k denotes the hydrodynamic stress tensor) 1982; Raupach, Coppin, & Legg, 1986) have associated this term with the work of drag on the interfacial boundary (term 17 in Eq. (17)), assuming the form-induced transport terms are insignificant. The exchange between the double-mean and form-induced flow fields φ Vm φ Tũiũj (∂ ū i /∂x j ) (term 5 in Eq. (19) and term 13 in Eq. (17)) is interpreted as mainly the conversion of the double-mean energy into to the form-induced energy (Yuan & Piomelli, 2014). The terms associated with the energy exchange at the boundary, i.e. the interfacial terms 14,15 in Eq. (19) and terms 11,12 in Eq. (21), are introduced in this paper and thus have never been assessed. Similarly, terms 3, 4, 5 in Eq. (17) and term 3 in Eq. (19), which appear in the equations due to potential correlations of the time-averaged quantities with the time porosity, are introduced here for the first time and their role remains to be clarified. At a first step, their assessment and interpretation can be done based on appropriate numerical simulations. The experimental capabilities will hopefully become available soon to complement simulations.
The equations proposed in this paper are applicable for flows over both mobile and immobile boundaries. For flows over immobile rough beds, the fluid averaging time T f is equal to the total averaging period T 0 , and thus φ T ≡ 1. The representation of the roughness geometry in the double-averaged equations then may be simplified by replacing φ VT and φ Vm with the roughness geometry function (or spatial porosity) φ, as φ VT = φ Vm = φ. The double-mean energy balance in Eq. (17) can therefore be rewritten as: Similarly, the form-induced energy balance in Eq. (19) takes the following form for immobile-bed flows: whereũ i = − ū i was used at the interface for modifying terms 14 and 15 in Eq. (19). The work of pressure on double-mean strain rate p ∂ ū i /∂x i (as ∂ ū i /∂x i is non-zero if φ spatially varies, Eq. (9)) and the work of pressure drag and viscous friction (surface integral terms) are terms that arise with opposite signs in both Eqs (25) and (26). Hence, these terms can be interpreted as bridges for energy exchange between the double-mean and form-induced fields due to the roughness spatial heterogeneity and the effect of drag. For fixed-bed flows, the energy from gravity is supplied to the double-mean flow only, while the energy to the form-induced filed is likely provided from the double-mean flow through work of the form-induced stresses on the double-mean strain rate − ũ iũj ∂ ū i /∂x j , as was mentioned above. The balance of the spatially-averaged turbulent energy for the fixed-bed flows follows from Eq. (21) as: As time porosity is φ T ≡ 1 and the mean strain rate is ∂ū i /∂x i = 0, the turbulent pressure-strain rate correlation (term 9 in Eq. (21)) vanishes to zero.

Application example: flow over an array of wind or tidal turbines
For some tasks, the equations proposed in Section 4 can be simplified by applying appropriate assumptions and boundary conditions. Such simplified equations may help clarify the flow interactions with mobile boundaries where energy is extracted from or supplied to the boundary motion. Below, we consider a hypothetical example related to flows over a large array of wind or tidal turbines, as represented schematically in Fig. 2. The characteristic length scale of the flow is assumed to be small enough for the Coriolis effects to be negligible. In such conditions, the flow is driven by a streamwise pressure gradient. The efficiency of the deployed array, in terms of power extraction from the flow, can be studied using the balances of the double-mean energy, form-induced energy, and spatiallyaveraged turbulent energy. For simplicity we use an averaging volume V 0 that includes all the turbines in the streamwise and spanwise direction, as noted in Fig. 2. In the vertical direction, Figure 2 Sketch of flow over a wind turbine array. The region within dotted lines represents the total averaging volume V 0 the averaging domain is thin enough to resolve hydrodynamic gradients in the surface-normal direction. We also assume flow two-dimensionality in terms of the double-averaged quantities, and so any spanwise derivatives of the double-averaged quantities are negligible, i.e. ∂ θ /∂x 2 = 0. Furthermore, we assume the following conditions: • The flow is steady, i.e. ∂ θ /∂t = 0, and fully developed uniform, i.e. ∂ θ /∂x 1 = 0. • Spatial correlations between time-porosity φ T and the timeaveraged fields are negligible, i.e. φ Tθ = 0.
• The flow is turbulent and fully developed and the effects of viscous fluid stresses on the double-averaged momentum balance are neglected, except viscous friction at the flow surface(s).
With these assumptions, the average loading f D from pressure drag and viscous friction on the turbine array can be expressed using the double-averaged momentum Eq. (15) as: where φ = φ VT is the space-time porosity. The equation for the double-mean energy Eq. (17) is simplified to the form: where viscous transport and viscous dissipation are assumed to be negligible, as in fully developed turbulent flows they are orders of magnitude smaller than their turbulent counterparts. The interpretation of terms in Eq. (29) is the same as in Eq. (17).
The equations for the form-induced and turbulent kinetic energy (Eqs (19) and (21)) can also be simplified, i.e.: Viscous transport effects are considered negligible and are not included in Eqs (30) and (31). Terms P D , P F and P T represent the work spent on the surface of the turbines by the doublemean, form-induced and turbulent velocity fields, respectively. The role of these terms is two-fold: (1) generation of the turbinescale turbulence (or "wake turbulence") by pressure drag, and (2) energy transfer with the rotating blades. If the turbine blades are "frozen" then the second role vanishes, as in the conventional situation of immobile rough-bed turbulent flows. Equations (28)-(31) may serve as a physically justified basis for numerical modelling of the flows through and over turbine arrays. Furthermore, through the interfacial terms P D , P F and P T they provide a way of coupling the turbulent flow with the turbine operations. It should be noted that equations, similar to Eq. (29), have already been used in studies of turbulent boundary layers over wind-turbine arrays to determine the mean power extracted by the turbines (Cal, Lebrón, Castillo, Kang, & Meneveau, 2010;Calaf, Meneveau, & Meyers, 2010). In these works, the rate of the energy extraction was fully assigned to term P D . It is worth emphasizing that in this respect Eq. (29) provides additional methodological support to Cal et al. (2010) and Calaf et al. (2010) who used equations originally developed for flows over immobile boundaries. Also, we believe that the terms P D are associated with actual energy extraction only partially, as in addition to them there may be contributions from the turbulent and form-induced terms P F and P T . However, no systematic information is available regarding turbine interactions with the turbulent and the form-induced fields in terms of the energy exchange. Thus, it is not clear how blade motions interact with the form-induced and turbulent flow fields. According to Eqs (29)-(31), the energy exchange between the turbine rotors and fluid motions (at the array scale) can occur from all three sources, i.e. the double-mean, forminduced and turbulent components of the flow. The energy exchange with these fields may be important when the total power that is extracted from the flow is determined. Considering similar boundary conditions for flows over tidal turbine arrays, Eqs (29)-(31) can be used for the study of such cases. Equations (29)-(31) can be applied to the analysis of data coming from high resolution experiments or simulations to investigate such problems.

Conclusions
We have proposed a set of the first-and second-order timeaveraged equations for mobile-boundary conditions that include mass, momentum, and momentum flux (i.e. stress) conservation equations. These equations can be used directly as a basis for consideration of mobile-boundary flows or in the spatially-averaged forms that are more suitable when dealing with flows over complex and mobile boundaries. Therefore, in addition to the time-averaged equations, we also proposed the double-averaged second-order equations that describe conservation of momentum fluxes and energy for the double-mean, form-induced, and turbulent components of the flow. The interrelations between these components as well as the effects of temporally and spatially varying roughness geometry are explicitly accounted for in the proposed equations. These features differentiate these equations from their counterparts that are already available for flows over immobile boundaries.
The complete interpretation of the terms of the derived balances of the double-mean, form-induced, and turbulent stresses and energies is not possible without involvement of highresolution data. As a first step, highly resolved data from direct numerical simulations of open-channel flows over mobile granular beds (Vowinckel & Fröhlich, 2012) were used in Vowinckel et al. (2017a) to assess the terms of the double-averaged momentum equation. The same dataset was then used in Papadopoulos et al. (2018a) to evaluate the double-mean and form-induced energy balances, and in Papadopoulos et al. (2018b) to assess the spatially-averaged turbulent energy balance.
Together with the double-averaged conservation equations for mass and momentum, the second-order double-averaged equations provide a solid framework for experimental and numerical investigations of hydraulic problems involving flow over or through mobile boundaries. This framework can help in significant data reduction and assist in understanding the key physical mechanisms involved in mobile-bed complex flows, at scales relevant to the selected averaging domains. The proposed energy budgets may provide useful insights into the interplay between boundary heterogeneity and mobility, flow (spatial) near-wall heterogeneity, and turbulence. After appropriate simplifications (as per the wind turbine example), the proposed equations can help in the conceptual development of improved hydraulic models and associated parameterizations. In addition to the design, optimization, and monitoring of wind turbine installations, applications include studying the effects of changing river bed morphology on the flow for flood control, the movement of fish schools, and river habitat management and restoration projects.  Figure A1 Sketch of a flow domain containing fixed particles (white) and rightward moving particles (shaded) at different time instants (a-e). The distribution of the indicator function γ (f) and its derivative ∂γ / ∂t with respect to time t at a position A is shown in (g) following Gray and Lee (1977).
(2) Then, in the same manner, time-averaging is performed on Eq. (6) to obtain the time-averaged continuity and momentum Eqs (7) and (8). (3) Equations (7) and (8) (5) Finally, the outcomes of steps 1 and 3 (based on Eqs (B1) and (B2)) are used to express the right-hand side time derivatives in Eq. (B3) to obtain the equation for local (at a point) turbulent stresses.

Double-averaged equation for form-induced stresses
(1) The starting point is to obtain the equation for the spatiallyaveraged mean stresses φ Vm φ Tūiūk by integrating Eq. (10) over V m , dividing by V 0 , and then applying the averaging theorems given in Eq. (A8). For brevity, the final equation is not shown in full here as it represents an intermediate step in the derivation.
(2) At a next step, the combined balance equation for the quantities φ VT ū i ū k , φ Vm φ Tũi ū k and φ Vm φ Tũk ū k is obtained by substituting the double-averaged mass and momentum conservation Eqs (14) and (15) in the identity: Eq. (15) (3) The third step involves the decompositionθ = θ +θ that is used to subdivide the spatially-averaged mean stress φ Vm φ Tūiūk into double-mean φ VT ū i ū k and forminduced φ Vm φ Tũiũk momentum fluxes as well as contributions from the porosity-velocity correlations, as explicitly shown in Eq.
(3). The balance equation for the form-induced