In this thesis we have developed a path towards large scale Finite Element simulations of turbulent incompressible flows. We have assessed the performance of residual-based variational multiscale (VMS) methods for the large eddy simulation (LES) of turbulent incompressible flows. We consider VMS models obtained by different subgrid scale approximations which include either static or dynamic subscales, linear or nonlinear multiscale splitting, and different choices of the subscale space. We show that VMS thought as an implicit LES model can be an alternative to the widely used physical-based models. This method is traditionally combined with equal-order velocity-pressure pairs, since it provides pressure stabilization. In this work, we also consider a different approach, based on inf-sup stable elements and convection-only stabilization. In order to do so, we define a symmetric projection stabilization of the convective term using an orthogonal subscale decomposition. The accuracy and efficiency of this method compared with residual-based algebraic subgrid scales and orthogonal subscales methods for equal-order interpolation is also assessed in this thesis. Furthermore, we propose Runge-Kutta time integration schemes for the incompressible Navier-Stokes equations with two salient properties. First, velocity and pressure computations are segregated at the time integration level, without the need to perform additional fractional step techniques that spoil high orders of accuracy. Second, the proposed methods keep the same order of accuracy for both velocities and pressures. Precisely, the symmetric projection stabilization approach is suitable for segregated Runge-Kutta time integration schemes. This combination, together with the use of block-preconditioning techniques, lead to elasticity-type and Laplacian-type problems that can be optimally preconditioned using the balancing domain decomposition by constraints preconditioners. The weak scalability of this formulation have been demonstrated in this document. Additionally, we also contemplate the weak imposition of the Dirichlet boundary conditions for wall-bounded turbulent flows. Four well known problems have been mainly considered for the numerical experiments: the decay of homogeneous isotropic turbulence, the Taylor-Green vortex problem, the turbulent flow in a channel and the turbulent flow around an airfoil.