The multiscale finite element method is applied in this paper to simulate groundwater flow in heterogeneous porous media. The method can efficiently capture the large-scale behavior of the solution without resolving all the small-scale features by constructing the multiscale finite element base functions that are adaptive to the local property of the differential operator. Both multiscale finite element method and conventional finite element method are applied to five 2-D and two 3-D groundwater flow problems, including a 2-D steady flow problem with continuous coefficients, a 2-D steady flow with gradual change in coefficients, a 2-D transient flow with gradual change in coefficients, a 2-D steady flow problem with an abrupt change in coefficients, a 2-D transient flow problem with an abrupt change in coefficients, and a 3-D steady and a 3-D transient flow problem with gradual change in coefficients in the horizontal direction and with abrupt change in the vertical direction, respectively. The applications demonstrate the main advantages of the multiscale finite element method, i.e., significantly reducing computational efforts and improving the accuracy of the solutions.