We present a residual-based turbulence model for problems with free surfaces. The method is derived based on variational multiscale ideas that assume a decomposition of the solution fields into overlapping scales that are termed as coarse and fine scales. The fine scales are further split hierarchically into fine-scales level-I and fine-scales level-II. The hierarchical variational problems that govern the two fine-scale components are modeled employing bubble functions approach. The model for level-II scales is variationally embedded in the mixed field level-I problem to yield a stable level-I formulation. Subsequently, the model for level-I scales that in fact constitutes the fine-scale turbulence model is then variationally injected in the coarse-scale variational form. A significant feature of the method is that it does not contain any embedded tunable parameters. To accommodate the moving boundaries we cast the formulation in an arbitrary Lagrangian–Eulerian frame of reference. The free surface boundary condition is imposed weakly which results in a formulation that conserves the volume of the fluid. A variety of benchmark problems show the accuracy and range of applicability of the proposed formulation and results are compared with published data. A wavy bed problem is investigated to show the interaction of turbulence generated at the bottom surface with the free surface thereby leading to irregular free surface elevations.