In this paper, we propose a semi-implicit approach for the time discretization of the Navier–Stokes equations with Variational Multiscale-Large Eddy Simulation turbulence modeling (VMS-LES). For the spatial approximation of the problem, we use the Finite Element method, while we employ the Backward Differentiation Formulas (BDF) for the time discretization. We treat the nonlinear terms arising in the variational formulation of the problem with a semi-implicit approach leading to a linear system associated to the fully discrete problem which needs to be assembled and solved only once at each discrete time instance. We solve this linear system by means of the GMRES method by employing a multigrid (ML) right preconditioner for the parallel setting. We validate the proposed fully discrete scheme towards the benchmark problem of the flow past a squared cylinder at high Reynolds number and we show the computational efficiency and scalability results of the solver in a High Performance Computing framework.