Modern measurement techniques have shown that property distributions in natural porous and fractured media appear highly irregular and nonstationary in a spatial statistical sense. This implies that direct statistical analyses of the property distributions are not appropriate, because the statistical measures developed will be dependent on position and therefore will be nonunique. An alternative, which has been explored to an increasing degree during the past 20 years, is to consider the class of functions known as nonstationary stochastic processes with spatially stationary increments. When such increment distributions are described by probability density functions (PDFs) of the Gaussian, Levy, or gamma class or PDFs that converge to one of these classes under additions, then one is also dealing with a so-called stochastic fractal, the mathematical theory of which was developed during the first half of the last century. The scaling property associated with such fractals is called self-affinity, which is more general that geometric self-similarity. Herein we review the application of Gaussian and Levy stochastic fractals and multifractals in subsurface hydrology, mainly to porosity, hydraulic conductivity, and fracture roughness, along with the characteristics of flow and transport in such fields. Included are the development and application of fractal and multifractal concepts; a review of the measurement techniques, such as the borehole flowmeter and gas minipermeameter, that are motivating the use of fractal-based theories; the idea of a spatial weighting function associated with a measuring instrument; how fractal fields are generated; and descriptions of the topography and aperture distributions of self-affine fractures. In a somewhat different vein the last part of the review deals with fractal- and fragmentation-based descriptions of fracture networks and the implications for transport in such networks. Broad conclusions include the implication that models based on increment distributions, while more realistic, are inherently less predictive than models based directly on stationary stochastic processes; that there is presently an unresolved ambiguity when a measurement is attempted in a medium that exhibits property variations on all scales; the strong possibility that log(property) increment distributions that appear to be described by the Levy PDF are actually superpositions of several PDFs of finite variance, one for each facies; that there are apparent similarities in the transport behavior of heterogeneous porous media and fractured rock at the field scale that appear to be related to the existence of a few preferential flow paths in both types of media; and finally, that additional carefully collected data sets are needed to clarify and advance the fractal-based theories, particularly in the case of three-dimensional fracture networks where few data are available. Further refinement is needed also in the understanding of instrument spatial weighting functions in heterogeneous media and how measurements in media exhibiting variations on all scales should be interpreted.