We address some unsolved methodological issues in modeling of natural rough-bed flows by critically examining existing approaches that parameterize rough-bed flows. These often utilize loosely defined variables. Here we suggest that double-averaging (in time and in a volume occupying a thin slab in a plane parallel to the mean bed) provides a rigorous, straightforward alternative that can aid in the parameterization process. We further present two examples: two-dimensional bed form and gravel bed flows. We argue that the double-averaging approach based on momentum equations should serve as a better methodological basis for modeling, phenomenological developments, and parameterizations. These equations explicitly include drag terms and form-induced momentum fluxes due to spatial heterogeneity of the time-averaged flow in the near-bed region. They also provide a solid basis for better definitions of basic flow variables including the shear stress partitioning into turbulent and form-induced momentum fluxes, skin friction, and pressure (form) drag. We show that for a range of rough-bed flows the vertical distribution of the double-averaged velocity consists of two distinct regions: (1) a linear region below roughness tops, and (2) a logarithmic region above them.