A new efficient multilevel upscaling procedure for single-phase saturated flow in porous media is presented. While traditional approaches to this problem have focused on the computation of an upscaled hydraulic conductivity, here the coarse-scale model is created explicitly from the fine-scale model through the application of operator-induced variational coarsening. This technique, which originated with robust multigrid solvers, has been shown to accurately capture the influence of fine-scale heterogeneous structure over the complete hierarchy of coarse-scale models that it generates. Moreover, implicit in this hierarchy is the construction of interpolation operators that provide a natural and complete multiscale basis for the fine-scale problem. Thus this new multilevel upscaling methodology is similar to the multiscale finite element method, and indeed, it attains similar accuracy in computations of the fine-scale hydraulic head and coarse-scale normal flux on a variety of problems; yet it is an order of magnitude faster on the examples considered here.