The full numerical simulation of turbulent flows in the context of industrial applications remains a very challenging problem. One of the difficulties is the presence of a large number of interacting scales of different orders of magnitude, ranging from the macroscopic scale to the Kolmogorov dissipation scale. In order to better understand and simulate these interactions, new algorithms of incremental type have recently been introduced, stemming from the theory of infinite-dimensional dynamical systems (see, e.g., the algorithms of Foias, Jolly et al. (1988), Foias, Manley and Temam (1988), Laminie et al. (1993, 1994), Marion and Temam (1989, 1990), and Temam (1990)). These algorithms are based on decompositions of the unknown functions into a large- and a small-scale component, one of the underlying ideas being to approximate the corresponding attractor. In the context of finite element methods, the decomposition of a solution into small- and large-scale components appears when we consider hierarchical bases (Yserentant, 1986; Zienkiewicz et al., 1982). In the present article we derive several estimates of the size of the structures for the linear and nonlinear terms which correspond to interactions of different hierarchical components of the velocity field, and also their time variation. The one-step time variation of these terms can be smaller than the expected accuracy of the computation. Using this remark, we implement an adaptive spatial and temporal multilevel algorithm which treats differently the small- and large-scale components of the flow. We derive several a priori estimates in order to study the perturbation introduced into the approximated equations. All the interaction terms between small- and large-scale structures are frozen during several time steps. Finally, we implement the multilevel method in order to simulate a bidimensional flow described by the Burgers equations. We perform a parametric study of the procedure and its efficiency. The gain in CPU time is significant and the trajectories computed by our multi-scale method remain close to the trajectories obtained with the classical Galerkin method.