Incompressible turbulent flow plays a key role in engineering and physics. It often requires the numerical resolution and modeling of multiple spatiotemporal scales. This work presents a numerical method that combines a residual-based variational multiscale (VMS) formulation with a half-explicit Runge-Kutta (RK) time integration method. The VMS framework can effectively address the spatial multiscale problem for the large eddy simulation of turbulent flows. Built upon this formulation, by using higher-order, high-continuity splines, we construct a scheme with favorable dissipation and dispersion properties [1]. The higher-order half-explicit RK methods ensure high temporal accuracy. For the incompressible Navier-Stokes equations, velocity and pressure can attain the theoretical order of accuracy [2]. Compared to fully implicit methods, it is computationally more efficient, as it only requires solving a linear Darcy-type problem at each stage. Moreover, in contrast to the explicit multistep methods, the RK methods offer better stability, have no start-up problem, and readily support adaptive time stepping. In particular, the subgrid-scale terms of the traditional residual-based VMS involve time derivatives. This renders the semi-discretized problem into an implicit system of ordinary differential equations, making it quite challenging to directly apply the RK method. To address this issue, we adopt the Rothe method [3] to discretize the problem in time first. The resulting Darcy-type problem is then spatially discretized using the VMS. Furthermore, we re-design the subgrid model for the Darcy-type problem within the framework of VMS. The proposed method is validated through two numerical benchmarks. First, a convergence study for the Taylor-Green vortex flow is performed. Second, numerical simulations of isotropic turbulence are investigated. The results confirm the viability of our method and demonstrate its high spatiotemporal accuracy.