Residual-based stabilized finite element (FE) techniques for the Navier–Stokes equations lead to numerical discretizations that provide convection stabilization as well as pressure stability without the need to satisfy an inf–sup condition. They can be motivated by using a variational multiscale (VMS) framework, based on the decomposition of the fluid velocity into a resolvable FE component plus a modelled subgrid-scale component. The subgrid closure acts as a large eddy simulation turbulence model, leading to accurate under-resolved simulations. However, even though VMS formulations are increasingly used in the applied FE community, their numerical analysis has been restricted to a priori estimates and convergence to smooth solutions only, via a priori error estimates. In this work, we prove that some versions of these methods (based on dynamic and orthogonal closures) also converge to weak (turbulent) solutions of the Navier–Stokes equations. These results are obtained by using compactness results in Bochner–Lebesgue spaces.