In this paper we present a numerical method, the finite volume multiscale finite element method (FVMSFEM), for solving the groundwater flow problems in heterogeneous porous media spanning over many scales. This method is based on an efficient coupling between the finite volume discretization and the multiscale finite element base functions. It can efficiently capture the large-scale structure of the solution on a coarse grid without resolving all the small-scale features and is locally conservative. The underlying idea is to estimate the macroscopic fluxes across the finite control volume interface segments, which bring the local small-scale information of the medium property to the large scales, by employing the multiscale finite element base functions. We describe the strategy for constructing such a method on the basis of the transient flow problems in porous media. Numerical experiments are carried out for groundwater flow in porous media with a random lognormal conductivity field to demonstrate the efficiency and accuracy of the developed method.