To extend the application of the variational multiscale method to ocean engineering, 2D residual-based large eddy simulations of flows around a circular cylinder close to a flat plane have been considered in the present study. A generalized-α method is applied to solving variational forms of the Navier–Stokes equations. Flows around a spanned marine pipeline at a subcritical Reynolds number 1.31×10⁴ have been computed. The hydrodynamic forces are compared to the published experimental data by a PIV and numerical results by a RANS and a classical LES. Although the magnitude of the predicted hydrodynamic forces are over-predicted compared to the experimental data, they are in good agreement with those obtained by the classical LES. The effect of the gap distance between the cylinder and the seabed on the hydrodynamic forces are analyzed. It is found that the drag coefficients and the root-mean-square of the lift coefficients decrease as the gap distance decreases. The flow characteristics around the spanned pipeline are examined. It is demonstrated that the proximity of the cylinder to the seabed severely affects the wake flow of the cylinder, and that the vortex shedding is suppressed as G∕D<0.3, which is in agreement with the experimental and other numerical results. The pressure distributions along the seabed appear to coincide with the experimental data. It is significant that the present method results in more reliable shear stresses (friction velocities) along the seabed than the RANS, which is very important for the research of the sediment transportation.