We study numerically and analytically the turbulent diffusion characteristics in a low-Froude-number turbulent shear flow beneath a free surface. In the numerical study, the Navier–Stokes equations are solved directly subject to viscous boundary conditions at the free surface. From an ensemble of such simulations, we find that a boundary layer develops at the free surface characterized by a fast reduction in the value of the eddy viscosity. As the free surface is approached, the magnitude of the mean shear initially increases over the boundary (outer) layer, reaches a maximum and then drops to zero inside a much thinner inner layer. To understand and model this behaviour, we derive an analytical similarity solution for the mean flow. This solution predicts well the shape and the time-scaling behaviour of the mean flow obtained in the direct simulations. The theoretical solution is then used to derive scaling relations for the thickness of the inner and outer layers. Based on this similarity solution, we propose a free-surface function model for large-eddy simulations of free-surface turbulence. This new model correctly accounts for the variations of the Smagorinsky coefficient over the free-surface boundary layer and is validated in both a priori and a posteriori tests.