Natural geological formations generally show multiscale structural and functional heterogeneity evolving over many orders of magnitude in space and time. In subsurface hydrological simulations the geological model focuses on the structural hierarchy of physical sub units and the flow model addresses the functional hierarchy of the flow process. Flow quantities like pressure, flux and dissipation relate to each other by constitutive relations and structural sub-unit parameters like porosity and hydraulic permeability. Hydraulic permeability includes the (steady state) intrinsic permeability of the solid phase and the (time dependent) relative permeability. The permeability is the dominant parameter and a highly heterogeneous parameter affecting the groundwater flow at field scale. Laboratory experiments provide measurements of the sub-unit parameters on a fine scale. If laboratory measurements are treated stochastically within the geological model, then the structural model of the subsurface should be built on the same scale as these indicating measurements. Fully resolved flow simulations on the field scale are however intractable and a new adaptive multiscale technique has therefore been developed. Though constitutive relations may change at different scales, Darcy’s law is supposed to remain valid on both laboratory scale and field scale. Nowadays upscaling methods are applied, which aim to propagate information over this hierarchy of scales both functionally and structurally from the fine scale to the coarse scale and not vise versa, by computing effective or equivalent material behavior. Permeability is not an additive parameter so it is not possible to calculate an equivalent coarse scale permeability as a simple average of fine scale measurements. A flow criterion or a criterion about energy dissipation often defines an equivalent permeability. Only if the scale of variation is much smaller than the coarse observation scale then the equivalent permeability matches the effective permeability. This effective permeability is a constant second order tensor variable on the coarse scale, whereas the equivalent permeability is non-unique and depends on the boundary conditions of the sample domain. The effective permeability holds for discrete hierarchical systems where scales can be decomposed. The newly developed adaptive multiscale technique extends the original two-level multiscale finite element method to a hierarchy of scales. Multiscale finite element methods capture the fine scale behavior on a coarse mesh by multiscale basis functions. The weights of the multiscale basis functions follow from local flow simulation. The procedure removes fine scale nodes from the subdomain, but introduces errors at the subdomain boundaries. The two-level method forms a class of subdomain decomposition techniques. It can be shown that a sequential implementation of the method is not faster than an optimal solver like the full multigrid solver, however the method is suitable for parallel implementation. The proposed method is based on a conformal nodal finite element formulation over simplex elements. The conformal finite element method obtains mass conservation on a nodal basis, and does not preserve continuity of flux over the inter-element boundaries. A mesh refinement criterion detects zones in which large errors occur over the element edges, and an adaptive refinement procedure enriches the mesh locally to correct the error in the velocity field. Multiscale basis functions follow from solving local flow problems over patches of simplex elements. Linear boundary conditions close oversampled subdomains and reduce the effect of the imposed boundary condition on the patch. This procedure obtains more accurate coarse scale behavior then a procedure that operates on the patch directly. However, over-sampled local flow problems introduce discontinuities in the basis functions and introduce new nodal connectivities on the coarser scale. For this reason closure of the local flow problems by dimensionally reduced flow problems is preferred. A second refinement criterion compares oversampled and non-oversampled function values. Pressure-dissipation averaging approximates the multiscale coarse grid operator and supports a functional adaptive formulation. The multiscale averaging procedure computes equivalent permeability tensor components, and reproduces the sparse matrix structure on the coarse scale. The loading cases for the local problems follow from a summation of multiscale basis functions. The multiscale basis functions extrapolate the coarse scale solution to the fine scale. On this scale discontinuities in the velocity field are detected and compared to the refinement criterion. The computed equivalent permeability is used in the framework of a geometric multigrid solver to compute the coarse scale operators. The hierarchy of multiscale basis functions, which relate the pressure on each coarse level to the next fine level, generates the intergrid transfer operators. The proposed approach provides a robust and efficient algorithm, based on the concept of the multiscale finite element method, for simulating partly saturated subsurface flow and fully coupled solute transport and heat transport through hierarchical heterogeneous formations. The multiscale finite element formulation produces numerically homogenized discrete flow equations, and upscales the permeability. The adaptive formulation obtains locally refined velocity fields, which support accurate transport computations. The performance of the method is illustrated by a set of realistic case studies.