Freestream turbulence in rivers is a key contributor to the flux of dissolved nutrients, carbon, and other ecologically important solutes into porewater. To advance understanding of turbulent hyporheic exchange and porewater transport, we investigate flow over and through a rough bed of spheres using large eddy simulation (LES). We apply double averaging (combined space and time averaging) to the LES results to determine the mean velocity distribution, momentum balance, and drag forces. Our simulations show large-scale freestream structures interacting strongly with vortices generated at the surfaces of individual spheres to control turbulent momentum fluxes into the bed. The transition between turbulent flow and Darcy flow occurs over the first row of spheres, where turbulence decays rapidly and turbulent kinetic energy, Reynolds stress, and drag forces peak. Below this region, turbulence is only present in the high-velocity flow in open pore throats. Experimental observations suggest that minimum mean porewater velocity occurs in the first open pore space below the transition region, but our results show that the minimum occurs between the first and second pore spaces. The simulation mean porewater velocities are approximately half those captured in measurements because the model resolves the entire flow continuum while measurements can access high-velocity fluid in open pores. The high-resolution dual time-space averaging of the LES resolves both turbulent and mean flow features that are important to interfacial solute and particle fluxes, providing a means to include turbulent hyporheic exchange in upscaled river models, which has not been achieved to date.