Numerical solution of unsaturated water flow in heterogeneous porous media is difficult not only because of the spatial variability in the parameters used to characterize the relevant physical properties of the natural porous media but also because of the mathematical problems that arise in dealing with the nonlinearities. In this study a multiscale finite element linearization scheme is presented for effectively simulating unsaturated flow in heterogeneous porous media spanning over many scales. The central goal of this method is to obtain the large-scale solution of Richards' equation with heterogeneous coefficients accurately and efficiently on a coarse grid without resolving all the small-scale details. The underlying idea is to employ the Slodicka linear relaxation approximation scheme to address the nonlinear characteristics of the equation and to use the multiscale finite element base functions to account for the spatial variability in the equation coefficients. We describe the principle for constructing such a method and give an algorithm for implementing it. Numerical experiments are carried out for the unsaturated flow equation with both periodic and randomly generated lognormal hydraulic conductivity to demonstrate the efficiency and accuracy of the proposed method.