Solute transport in the subsurface is strongly affected by the heterogeneity of the porous medium at all scales. In this context, (macro)dispersion denotes the process of a solute cloud expanding its spatial extent upon transport (e.g., Dagan et al., 2003;Gelhar et al., 1979). In heterogeneous porous media, there are two contributions to dispersion: advective spreading, which is the volume-preserving distortion of the plume due to the non-uniform flow field, and diffusive mixing, which is the dilution of the plume by Brownian motion (Dentz et al., 2011).Spreading is deterministically caused by the geometry of the porous medium and the distribution of hydraulic conductivity. It leads to an increase of the plume surface while the plume volume remains constant. By contrast, mixing increases the volume over which the plume is distributed, implying a reduction of concentration in the plume center (Kapoor & Kitanidis, 1996). With the increased surface area of the plume and steepening of concentration gradients due to contraction of fluid elements perpendicular to the direction of their stretching (Dentz et al., 2011;Le Borgne et al., 2015), advective spreading in heterogeneous media fosters enhanced diffusive mixing. However, over very long periods of time, spreading increases much faster than mixing, so that applying parameterizations that are good for advective spreading to solute mixing leads to a strong overestimation of mixing.Only mixing facilitates chemical reactions of reactants that are initially separated. Therefore it is of vital interest to separate the contributions of spreading and mixing to overall dispersion. Mixing is an irreversible process, whereas the deformation of a solute plume due to spatially variable advection is perfectly reversible. In thermodynamics, irreversibility implies an increase of entropy. Along these lines, Kitanidis (1994) introduced the dilution index as metric of solute mixing. The dilution index is defined as the exponent of the Shannon entropy of the spatial concentration distribution. With the exception of very simple cases, such as a truly Gaussian plume in a homogeneous flow field, the dilution index cannot be computed analytically. In heterogeneous aquifers with

The adequate parameterization of solute transport in natural, inherently heterogeneous formations has been the subject of intensive research in hydrogeology and environmental engineering over the last four decades. Classical Fickian macrodispersion concepts, parameterizing the effects of unresolved spatial variability on field-scale transport by a diffusion-type macrodispersion term with constant coefficients, can neither reproduce the evolution of solute spreading over shorter distances and travel times nor can it predict long tailing observed in field investigations (e.g., Adams & Gelhar, 1992;Haggerty et al., 2001;Kang et al., 2015). First-order perturbative approaches, typically assuming second-order stationary multi-Gaussian log-hydraulic conductivity (ln K) fields and uniform-in-the-mean hydraulic gradients, are good in predicting second-central spatial moments of solute plumes as long as the variance of ln K remains small (