In this paper we investigate the large-eddy simulation (LES) of the interaction between a turbulent shear flow and a free surface at low Froude numbers. The benchmark flow field is first solved by using direct numerical simulations (DNS) of the Navier–Stokes equations at fine (1282 × 192 grid) resolution, while the LES is performed at coarse resolution. Analysis of the ensemble of 25 DNS datasets shows that the amount of energy transferred from the grid scales to the subgrid scales (SGS) reduces significantly as the free surface is approached. This is a result of energy backscatter associated with the fluid vertical motions. Conditional averaging reveals that the energy backscatter occurs at the splat regions of coherent hairpin vortex structures as they connect to the free surface. The free-surface region is highly anisotropic at all length scales while the energy backscatter is carried out by the horizontal components of the SGS stress only. The physical insights obtained here are essential to the efficacious SGS modelling of LES for free-surface turbulence. In the LES, the SGS contribution to the Dirichlet pressure free-surface boundary condition is modelled with a dynamic form of the Yoshizawa (1986) expression, while the SGS flux that appears in the kinematic boundary condition is modelled by a dynamic scale-similarity model. For the SGS stress, we first examine the existing dynamic Smagorinsky model (DSM), which is found to capture the free-surface turbulence structure only roughly. Based on the special physics of free-surface turbulence, we propose two new SGS models: a dynamic free-surface function model (DFFM) and a dynamic anisotropic selective model (DASM). The DFFM correctly represents the reduction of the Smagorinsky coefficient near the surface and is found to capture the surface layer more accurately. The DASM takes into account both the anisotropy nature of free-surface turbulence and the dependence of energy backscatter on specific coherent vorticity mechanisms, and is found to produce substantially better surface signature statistics. Finally, we show that the combination of the new DFFM and DASM with a dynamic scale-similarity model further improves the results.