We discuss the development of a two-level iterative solver for Stokes flow simulations in porous media, focusing on its application to compute absolute permeability from pixel-based microstructures. The Schur complement plays a key role in converting the saddle-point system arising from the weak form of the Stokes equations into a condensed system for velocity, solvable using a robust and efficient two-level solver. At the first (external) level, we employ the conjugate gradient method, though this choice is not restrictive, as alternatives like the minimal residual method can also be accommodated. As in many short-recurrence Krylov subspace methods, the global matrix of the condensed system does not need to be stored in memory, as long as a function to compute matrix-vector products is available. The inverse matrix from the Schur complement is replaced by subsystem solutions during the implementation of matrix-vector products, similarly to some preconditioning techniques. These subsystems define the second (internal) level, which can also be handled iteratively. In addition, the proposed solver is designed within the framework of image-based simulations via the finite element method, making it suitable for memory-efficient, matrix-free implementations on structured meshes. We demonstrate the solver's applicability through numerical results on a range of image-based porous media characterization problems, observing stable and fast convergence patterns in terms of iterations. Among these, we specifically apply the solver to characterize two distinct porous media imaged via X-ray micro-computed tomography: a fiber-reinforced composite and a sandstone sample.