We develop a residual-based variational multiscale (RBVMS) method based on isogeometric analysis for large-eddy simulation (LES) of wind-driven shear flow with Langmuir circulation (LC). Isogeometric analysis refers to our use of NURBS (Non-Uniform Rational B-splines) basis functions which have been proven to be highly accurate in LES of turbulent flows (Bazilevs, Y., et al. 2007, Comput. Methods Appl. Mech. Eng., 197, pp. 173–201). LC consists of stream-wise vortices in the direction of the wind acting as a secondary flow structure to the primary, mean component of the flow driven by the wind. LC results from surface wave-current interaction and often occurs within the upper ocean mixed layer over deep water and in coastal shelf regions under wind speeds greater than 3 m s−1. Our LES of wind-driven shallow water flow with LC is representative of a coastal shelf flow where LC extends to the bottom and interacts with the sea bed boundary layer. The governing LES equations are the Craik-Leivobich equations (Tejada-Martínez, A. E., and Grosch, C. E., 2007, J. Fluid Mech., 576, pp. 63–108; Gargett, A. E., 2004, Science, 306, pp. 1925–1928), consisting of the time-filtered Navier-Stokes equations. These equations possess the same structure as the Navier-Stokes equations with an extra vortex force term accounting for wave-current interaction giving rise to LC. The RBVMS method with quadratic NURBS is shown to possess good convergence characteristics in wind-driven flow with LC. Furthermore, the method yields LC structures in good agreement with those computed with the spectral method in (Thorpe, S. A., 2004, Annu. Rev. Fluids Mech., 36, pp. 584 55–79) and measured during field observations in (D’Alessio, S. J., et al., 1998, J. Phys. Oceanogr., 28, pp. 1624–1641; Kantha, L., and Clayson, C. A., 2004, Ocean Modelling, 6, pp. 101–124).