Random walks are a natural starting point for modeling solute migration through fractured rock. In special cases, they lead to the common advection-dispersion equation (ADE), which models a multiGaussian limit. In a more general setting, the particles may be restricted to moving with a small number of preferred directions. The fractures may also lead to non-Fickian scaling that is different in any of the preferred fracture directions. Particles may also "stick" in the same spot for random, heavy-tailed amounts of time. These apparent complications are easily handled by the same mathematical machinery that gave us the ADE. The spatial dispersion derivative may contain a fractional-order Laplacian and the time derivative may also be of fractional order less than unity. The order of the Laplacian is actually a matrix that describes the scaling rate in the fracture directions. To compare, the traditional ADE is of order 2I, where I is the identity matrix, and the scaling rate (the Hurst matrix) is the inverse of this.

The key features of the model are the ease of analytic solution, the heavy tails of the breakthrough curves, and the intuitive shapes (boomerangs, rabbit ears) that the plumes may attain due to the influence of a small number of fracture sets. The model is parsimonious because the geometry and scaling of the media are part of the governing equation. We show examples of this behavior in fractured and granular aquifers.