Fluid flow in porous media is a multiscale process where the effective dynamics, which is often the goal of a computation, depends strongly on the porous micro structure. Resolving the micro structure in the whole porous medium can, however, be prohibitive. Novel numerical methods that efficiently approximate the effective flow but resolve only a carefully selected reduced portion of the porous structure are of great interest. In this thesis we propose new numerical multiscale methods for Stokes flow in two- and three-scale porous media. First, we propose the Darcy--Stokes finite element heterogeneous multiscale method (DS-FE-HMM). The method is based on solving the Darcy equation on a macroscopic mesh using the finite element method with numerical quadrature, where the unknown permeability is recovered from micro finite element solutions of Stokes problems that are defined in sampling domains centered at macroscopic quadrature points. An adaptive scheme based on a posteriori error analysis is proposed, where micro-macro mesh refinement is driven by residual-based indicators that quantify both the micro and macro errors. Second, to address the increasing cost of solving the micro problems as the macroscopic mesh is refined, we combine the DS-FE-HMM with reduced basis (RB) method and propose a new multiscale method called the RB-DS-FE-HMM. Efficiency and accuracy of the method relies on a parametrization of the micro geometries and on the Petrov-Galerkin RB formulation that provides a stable and fast evaluation of the effective permeability. A residual-based adaptive mesh refinement scheme is proposed for the macroscopic problem. To achieve a conservative approximation we also combine and analyze a coupling of the RB method with a different macroscopic scheme based on the discontinuous Galerkin finite element method (DG-FEM). Finally, we consider a three-scale porous media model with macro, meso, and micro scale. At the intermediate meso scale the medium is composed of fluid and porous parts and the fluid flow is modeled with the Stokes-Brinkman equation. A three-scale numerical method is derived and an efficient algorithm based on the RB method and empirical interpolation method on the micro and meso scale is proposed.