The variational multiscale method provides a methodical framework for large eddy simulation of turbulent flows. In this work, a particular implementation in the form of a three-level finite element method separating large resolved, small resolved, and unresolved scales is proposed. Residual-free bubbles are used for the numerical approximation of the small-scale momentum equation. A stabilizing term is added, in order to take into account the effect of the small-scale continuity equation. This implementation guarantees the stability of the method without further provisions and offers substantial computational savings on the small-scale level. Furthermore, it is accounted for the unresolved scales by a specific dynamic modelling procedure. The method is tested for two different turbulent flow situations. Copyright © 2005 John Wiley & Sons, Ltd.