This paper introduces an hp-refinement high-order finite element discretization based on a modified hierarchical Jacobian polynomial basis into a residual-based large eddy simulation (LES) by variational multiscale method (VMM). A high-order finite element basis is achieved by means of p-refinement by introducing a modified hierarchical Jacobian polynomial basis. The benchmark turbulent channel flow is computed using in-house finite element codes with three different orders of the Jacobian polynomial basis. The near-wall flow field and turbulent statistical results obtained using the same degrees of freedom (DOFs) are compared. The results show that the accuracy is significantly improved with the increasing order. The transportation equation of the turbulent kinetic energy is then derived with the VMM, and statistical dissipations are presented to clarify the improvement in the numerical accuracy by high-order simulations. In addition, the matrix structures of linear systems built using different orders of polynomial basis are presented, and the variation in the computational cost with the order of the finite element basis is explained. Furthermore, comparison of the computational costs for cases with equivalently accurate results is made to assess the computational efficiency of the high-order numerical method. Consequently, the results provide compelling evidence for the power of the high-order LES technique with the hp-refinement discretization, which not only reduces the DOFs to be solved, but also achieves the higher computational efficiency than low-order methods.