1 Introduction

Transport of a conservative solute through spatially heterogeneous geological formations is an important and challenging research area in hydrogeology and environmental engineering. Typically modeled by virtue of the Advection–Dispersion Equation (ADE), a major challenge is the representation of these heterogeneous formations. Since a full representation is made impossible by both numerical and measurement limitations, effective descriptions are used which are derived at the practically relevant scale of interest (Renard and De Marsily 1997; Cushman et al. 2002; Lunati et al. 2002; Frippiat and Holeyman 2008).

The traditional framework assumed in such studies is that the scales of heterogeneity can be clearly separated. These scales can range from the pore scale to the quasi-homogeneous meso or Darcy scale and eventually to the macro scale (Bear 1972; Dagan 1989; Berkowitz 2002; Neuman 2005). This assumption of clearly separated scales facilitates the definition of a representative-elementary volume (REV), which in turn allows the derivation of the coefficients of the ADE on the relevant scales (Whitaker 1999).

Within this picture of clearly separated scales, diffusion is the relevant process at the pore scale, whereas at the quasi-homogeneous meso scale, solute spreading is modeled by dispersion. This is necessary to account for the effects of pore-scale velocity variations, otherwise lost by the upscaling process (Scheidegger 1954; Bear 1972). To parameterize the effects of the heterogeneities at the macro scale, it is common practice to apply macro-scale or macrodispersion tensors (Sudicky 1986a; Woodbury and Sudicky 1991; Boggs et al. 1992; Gelhar et al. 1992). Compared to the lab-scale or meso-scale dispersion tensor, the macrodispersion tensor, or at least its longitudinal component, is typically several orders of magnitude larger due to comprising additionally the effects of macro-scale velocity variations. This experimentally described macro-scale property is modeled mathematically by two different approaches: the ensemble as well as the effective dispersion (Cirpka and Kitanidis 2000; Dentz et al. 2011).

The ensemble dispersion tensor \({\mathbf{D}}^{\mathrm{ens}}\) is defined as half the rate of change of the second central moments of the ensemble-averaged concentration distribution. It describes the dispersion properties of a (hypothetical) ensemble of heterogeneous aquifers as a whole, including the uncertainty in the center of mass. The effective dispersion tensor \({\mathbf{D}}^{\mathrm{eff}} \) is defined as the ensemble-average of half the rate of change of the second central moments of the spatial concentration distributions in a single realization of a heterogeneous aquifer. This quantity characterizes the expected plume spread in a typical realization of the medium. In general, the ensemble dispersion tensor is easier to calculate but, due to the aforementioned additional uncertainty, the effective dispersion tensor is for most applications a better representation of the experimentally observed dispersion. Since it is often assumed that the difference between the two can be neglected, the ensemble dispersion tensor is often used (Kabala and Sposito 1991; Shvidler 1993; Fiori 1996).

Yet for heterogeneous media with separated scales, the two dispersion tensors are time-dependent as well as different for a transient period during which the solute transport is called anomalous or non-Fickian (Attinger et al. 1999; Dentz et al. 2000a, b, 2002). This transient behavior, however, vanishes at larger times and the system is said to be ergodic. The length of this pre-asymptotic behavior depends both on the type (ensemble or effective dispersion) as well as on the correlation length of the heterogeneities. The aforementioned practice of neglecting the differences between ensemble and effective dispersion is therefore error prone during this early period of time, but becomes justified eventually.

But this aforementioned picture of a subsurface with heterogeneities on a discrete hierarchy of clearly separated scales has been questioned by numerous studies that advocate for a continuous hierarchy of scales of heterogeneities instead (Pickens and Grisak 1981; Neuman 1990; Kemblowski and Wen 1993; Desbarats and Bachu 1994; Fiori 1998; Kim et al. 2004). Since the heterogeneities appear with similar characteristics on different scales such media have the properties of fractals and have consequently been called fractal media (Wheatcraft and Tyler 1988; Sahimi 1993; Zhan and Wheatcraft 1996; Bonnet et al. 2001; Molz et al. 2004).

Such a conceptualization of the subsurface as a fractal medium poses several questions with respect to the results derived under the assumption of separated scales, i.e. for a non-fractal medium. In non-fractal media, the ensemble and effective dispersion coefficients remain time-dependent and different for a transient period, the length of which depends on the correlation length of the heterogeneities. Since fractal media do not have a finite correlation length, one could expect this non-Fickian behavior to prevail indefinitely. Past studies tentatively confirm this notion. For example, both Neuman (1995) and Bellin et al. (1996) could show such indefinite non-Fickian behavior for the longitudinal dispersion, while neglecting local dispersion. In a subsequent study Fiori (2001) could derive semi-analytical solutions for the asymptotic time-dependent behavior of the longitudinal dispersion in two-dimensional media accounting for local dispersion. Depending on the fractal regimes Fiori (2001) could derive different time-dependent behavior showing anomalous behavior at infinite times. Attinger et al. (2004), neglecting again local dispersion, showed for a two-dimensional medium, that the transverse ensemble and effective dispersion coefficients approach zero, due to the incompressibility of the velocity field. In contrast to that, the results presented by Di Federico and Neuman (1998a) revealed little difference between the temporal behavior of the dispersion coefficient in media with a truncated fractal compared to media with a classical exponential variogram function. This surprising similarity between fractal and classical media was due to the aforementioned truncation of the fractal behavior at larger scales, which made their variograms (see Di Federico and Neuman 1998b) as well as the resulting transport behavior very alike. More recently Pannone (2017) investigated the impact of the Péclet number of the longitudinal macrodispersion coefficient, showing that both Fickian and non-Fickian behavior can be found depending on the fractality of the medium. Despite this large number of studies, there is no study to this date, where longitudinal and transverse dispersion are jointly investigated in a fully three-dimensional medium and explicit descriptions for the dispersion coefficients are presented. The gap of knowledge is especially relevant with respect to the transverse dispersion coefficient, since its importance has been emphasized by numerous studies (Cirpka et al. 1999, 2011; Klenk and Grathwohl 2002; Olsson and Grathwohl 2007; Chiogna et al. 2011; Rolle et al. 2012).

A second and perhaps more general question is the type of the macro-scale transport model, i.e., the validity of the ADE itself (Sposito et al. 1986; Neuman and Tartakovsky 2009). The basic requirement for the use of such an effective single-continuum model is the definition of an REV, wherein complete small-scale mixing of the relevant quantities is assumed. However, this requirement is not met in fractal media, since the different scales interact seamlessly. Modeling schemes that account for this problem are for example Multi-Rate Mass Transfer (MRMT) models (Haggerty and Gorelick 1995) or the Continuous-Time-Random-Walk (CTRW) (Berkowitz et al. 2006). Within the context of single-continuum models, trying to amend the classical ADE-based approach, two different approaches can be distinguished. In the first approach the concept of the ADE remains unchanged but the incomplete small-scale mixing is represented by time-dependent parameters (Dagan 1988; Fiori 2001; Sanchez-Vila et al. 2010). An alternative approach is the use of a fractional ADE (fADE), which is able to describe anomalous non-Fickian transport by virtue of fractional derivatives in space or time (Huang et al. 2006; Zhang et al. 2007; Gao et al. 2009; Schumer et al. 2009). Whereas in case of separated scales the anomalous behavior is limited to a transient period, this behavior might be found in fractal media at all times (Fiori 2001). It therefore stands to reason to prefer a fADE in case of solute transport in fractal media.

To investigate these questions, the remainder of this paper is organized as follows: in Sect. 2, we describe our model and the applied methods for the investigation of the dispersion coefficients in fractal media. In Sect. 3, we present the explicit expressions for the ensemble and effective dispersion coefficients in fractal media and the results for the temporal behavior of these quantities. The paper ends with conclusions in Sect. 5.

2 Model and methods

In this study we use the following notations: vector and tensor quantities are represented by bold characters, all quantities related to fractal media are marked by a hat \(\widehat{\left[ \cdot \right] }\) in order to distinguish them from their non-fractal counterparts, the overbar \(\overline{\left[ \cdot \right] }\) is used for the mean of heterogeneous quantities or an ensemble of realizations and a tilde \(\tilde{\left[ {\cdot } \right] }\) is used to denote the fluctuations about this mean value.

2.1 Transport equation

The time evolution of solute transport is typically modeled by the ADE (e.g. Dagan 1989; Gelhar 1993)

$$\begin{aligned} \frac{\partial }{\partial t} c( {\mathbf {x}}, t) + {\mathbf {u (x)}}\cdot \nabla c( {\mathbf {x}}, t) - {\mathbf { \nabla }} {\mathbf {D}} \nabla c({\mathbf {x}}, t) = \rho ({\mathbf {x}}) \delta (t). \end{aligned}$$
(1)

Here, \(c( {\mathbf {x}}, t)\) is the concentration of the transported solute, \({\mathbf {u (x)}}\) is the incompressible velocity and \({\mathbf {D}}\) is the local dispersion tensor. In the following analysis, we assume a point injection, implying \(\rho ({\mathbf {x}}) = \delta ({\mathbf {x}})\). As auxiliary conditions we assume an at least exponentially vanishing concentration at infinity for the unbounded domain.

The first coefficient in Eq. (1) is the velocity field \({\mathbf u(x)}\). In saturated heterogeneous media it can be computed by the Darcy equation

$$\begin{aligned} {\mathbf u (x) } = - \frac{K ({\mathbf {x}})}{\phi } \nabla h ({\mathbf {x}}). \end{aligned}$$
(2)

Here, \(h ({\mathbf {x}})\) is the hydraulic head and \(\phi \) is the porosity. Within the stochastic framework used in this study, the heterogeneous hydraulic conductivity \(K ({\mathbf {x}}) \) is assumed to be a translational invariant random field, i.e. the field is statistically homogeneous. Geostatistical investigations suggest K to be log-normally distributed (Freeze 1975; Sudicky 1986b). As a result, the log-conductivity \(Y = \ln (K)\) is normally distributed, i.e. it is a Gaussian random field.

The second coefficient in Eq. (1) is the local dispersion tensor \({\mathbf {D}}\) (Scheidegger 1961). Within the context of this study, we assumed this quantity to be spatially uniform.

2.2 Stochastic approach

The velocity field in Eq. (1) is identical to the center-of-mass velocity, i.e. the rate of change of the first moment of the moving solute cloud. By contrast, the local dispersion tensor \({\mathbf {D}}\) in Eq. (1) and the phenomenologically observable dispersion tensor \({\mathbf {D}}^{\mathrm{obs}}\) differ. The latter quantity can be estimated from the temporal development of the second central moments of the solute cloud:

$$\begin{aligned} D_{jn}^{\mathrm{obs}}(t) = \frac{1}{2} \frac{\mathrm{d}}{\mathrm{d}t} \left\{ m_{jn}^ {(2)}(t) - m_{j}^{(1)} (t) \;\, m_{n}^{(1)} (t) \right\} . \end{aligned}$$
(3)

Here \(m_{j}^{(1)} (t)\) and \(m_{jn}^{(2)} (t)\) are the first two raw moments of the properly normalized spatial concentration distribution in \(d\)-dimensions:

$$\begin{aligned} m_j^{(1)} (t)= & {} \int d {\mathbf {x}}\,\, x_j c ( {\mathbf {x}}, t), \end{aligned}$$
(4)
$$\begin{aligned} m_{jn}^{(2)} (t)= & {} \int d {\mathbf {x}}\,\, x_j x_n c ( {\mathbf {x}}, t). \end{aligned}$$
(5)

Using the stochastic approach means that the heterogeneous conductivity field \(K ({\mathbf {x}})\) of a single aquifer is a spatial random field. By virtue of Eqs. (1), (2) and (3), the concentration field \(c( {\mathbf {x}}, t)\), the velocity field \({\mathbf u(x)}\) and the associated observed dispersion \({\mathbf {D}}^{\mathrm{obs}}(t)\) become random spatial fields, too.

As outlined in the introduction, there are two ways to estimate the observable dispersion from an ensemble of realizations of \(c( {\mathbf {x}}, t)\). First, the so called ensemble dispersion coefficient

$$\begin{aligned} D_{jn}^{\mathrm{ens}}(t) = \frac{1}{2} \frac{\mathrm{d}}{\mathrm{d}t} \left\{ \overline{m_{jn}^{(2)}(t)} - \overline{m_{j}^{(1)}(t)} \ \overline{m_{n}^{(1)}(t)} \right\} = D_{jn} + \mathscr {D}_{jn}^{\mathrm{ens}}(t), \end{aligned}$$
(6)

which is defined by the second central moments of the averaged concentration distribution. Second, the effective dispersion coefficient

$$\begin{aligned} D_{jn}^{\mathrm{eff}}(t) = \frac{1}{2} \frac{\mathrm{d}}{\mathrm{d}t} \left\{ \, \overline{m_{jn}^{(2)}(t) - m_{j}^{(1)} (t) m_{n}^{(1)} (t)} \right\} = D_{jn} + \mathscr {D}_{jn}^{\mathrm{eff}}(t), \end{aligned}$$
(7)

which defined by the ensemble-average of the central second moments of the spatial concentration distributions. As can be seen from Eqs. (6) and (7), these estimates are conceptualized as a superposition of the local dispersion tensor \({\mathbf {D}}\) and some additional contributions \({\mathbf {\mathscr {D}}}^{[\cdot ]}\). This approach is therefore analogous to the notion of local dispersion, which is conceptualized as a superposition of molecular diffusion and mechanical dispersion. Like mechanical dispersion, which is caused by the heterogeneous flow paths at the pore scale, the contributions of \({\mathbf {\mathscr {D}}}^{[\cdot ]}\) are caused by the heterogeneous flow paths at the continuum scale. Both lead to an increase in the spreading of the solute at the scale of observation and have therefore to be considered.

2.3 Correlation function for isotropic fractal media

In case of non-fractal media, Gaussian correlation functions of the log-hydraulic conductivity field Y are often used, due to their simple mathematical treatment (Dagan 1988; Gelhar 1993; Dentz et al. 2000a; Zech et al. 2012, 2016). In a d-dimensional isotropic medium they have the form

$$\begin{aligned} C^{\mathrm{Y}} ({\mathbf {x}}) = \sigma _{\mathrm{Y}}^2 \exp \left( - \frac{1}{2} \frac{|{\mathbf {x}}|^2}{l^2} \right) , \end{aligned}$$
(8)

with l indicating the correlation length scale and \(\sigma _\mathrm{Y}^2 \) the variance of \( Y({\mathbf {x}}) \) (Dentz et al. 2000a). In fractal media, having no definite correlation length, algebraic functions like power-law functions are used

$$\begin{aligned} \widehat{C}^{\mathrm{Y}}({\mathbf {x}}) = \sigma _\mathrm{Y}^{2}z^{-\beta }({\mathbf {x}}), 0\le \beta \le 2. \end{aligned}$$
(9)

It should be noted that the integral length I, defined as the integral over the correlation function (Dagan 1986), is finite for Eq. (9) if \(1 < \beta \le 2\) holds. As a result, we can discern two different regimes: (1) for \(0 \le \beta \le 1\) where we have a true fractal regime without a finite integral length and (2) for \(1 < \beta \le 2\), where we have a quasi fractal regime, with a power-law correlation function but a finite integral length. The different transport behavior of these two regimes has been noted before (Dagan 1994; Bellin et al. 1996; Fiori 2001) and will be highlighted below, when discussing the results. The variable z in Eq. (9) is related to the Euclidean distance according to

$$\begin{aligned} z({\mathbf {x}}) = \sqrt{ \left( 1 + \frac{|{\mathbf {x}}|^{2}}{\widehat{l}^2}\right) }. \end{aligned}$$
(10)

The length scale \(\widehat{l}\) in Eq. (10) changes the behavior of the correlation function in Eq. (9) such that the singularity at \(|{\mathbf {x}}|=0\) is avoided, whereas the behavior for large distances with \(|{\mathbf {x}}| / l \gg 1\) is similar to that of a purely fractal medium as used by e.g. Dagan (1994), Bellin et al. (1996), Fiori (2001) and Pannone (2017). Real subsurface formations are always connected to some physical length scales, be it the total size of the aquifer or the support scale induced by the measurement process or the modeling assumptions (Neuman 1994; Neuman and Tartakovsky 2009). Consequently, this mathematical model is arguably a better physical representation than a pure fractal and has consequently found applications in the literature (Dagan 1994; Suciu 2010; Suciu et al. 2015).

The exponent \(\beta \) in Eq. (9) is related to the Hurst coefficient H and by extension to the fractal dimension \(d_f\) of the medium

$$\begin{aligned} \beta = d - d_f + 1 = d + H - 1. \end{aligned}$$
(11)

To make use of previous results by Attinger et al. (1999) and Dentz et al. (2000a), we use in this study a representation of the fractal correlation function, constructed as a superposition of Gaussian functions. According to Gradshteyn and Ryzhik (2007) (p. 370, expression 3.478 (1.)) we can reformulate Eq. (9) as follows

$$\begin{aligned} \widehat{C}^{\mathrm{Y}} ({\mathbf {x}}) = \sigma _\mathrm{Y}^2\frac{2}{\varGamma (\frac{\beta }{2})} \int \limits _0^{\infty } d\lambda \ \lambda ^{\beta -1} \ \exp \left[ -\left( 1 + \frac{ |{\mathbf {x}}|^{2} }{l^2}\right) \lambda ^2\right] . \end{aligned}$$
(12)

Using the following short hand notations

$$\begin{aligned} \int \limits _0^{\infty }d\varLambda = \frac{2}{\varGamma (\frac{\beta }{2})} \int \limits _0^{\infty } d \lambda \ \lambda ^{\beta -1}\ \exp (- \lambda ^2) \end{aligned}$$
(13)

and

$$\begin{aligned} \widehat{L}= \frac{l}{ \sqrt{2}\ \lambda }, \end{aligned}$$
(14)

we finally get

$$\begin{aligned} \widehat{C}^{\mathrm{Y}} ({\mathbf {x}}) = \sigma _{\mathrm{Y}}^2\int \limits _0^{\infty } d \varLambda \ \exp \left( -\frac{1}{2}\frac{ |{\mathbf {x}}|^{2} }{\widehat{L}^2}\right) . \end{aligned}$$
(15)

Due to the structural similarities between Eq. (8) and the integrand in Eq. (15), we can transfer the results obtained for non-fractal fields to the case of fractal fields, by describing the resulting fractal dispersion coefficients as a superposition of their non-fractal counterparts derived by Dentz et al. (2000a). Please note that by following Dentz et al. (2000a), we are using a first-order analysis of the resulting solute behavior. We discuss the ramifications of this, as well as other assumptions, at the end of the manuscript.

3 Results

In the following we present and discuss the results for the transport coefficients in isotropic fractal media with a point-like injection. The coefficients were evaluated for both the ensemble and effective quantities according to Attinger et al. (1999) and Dentz et al. (2000a) as introduced above, and both quantities were distinguished by their longitudinal and transverse components as explained in the introduction. For simplicity, we restrict the analysis to the case of an isotropic local dispersion tensor, \({\mathbf {D}} = {\mathbf {I}}D\), in which D is the local isotropic dispersion coefficient, and \({\mathbf{I}}\) is the identity matrix.

First, we give the full explicit analytical expressions for these four coefficients. These expressions will form the basis for the following analysis and discussion.

Next, we present results where we investigated the four different expressions numerically. Since the equations have no closed-form solutions, their evaluations were performed using the software package MAPLE. All expressions were evaluated with a mean groundwater velocity of \( \bar{u} = 1 m / d \), a flow factor of \( \gamma = \exp \left( \sigma ^2_Y/6 \right) \) and a variance of \(\sigma ^2_Y = 1\) for an adequate description of the effective conductivity and mean velocity for a mildly heterogeneous medium (Gelhar 1993). We investigated the response of the dispersion coefficients to variations in different physical variables such as the local dispersion coefficient D, the degree of fractality \(\beta \) and the characteristic length \(\widehat{l}\). In the graphical representation of these numerical results we use in the following a convention where the ensemble quantities are indicated by dashed lines, whereas the effective coefficients are denoted by solid ones.

After presenting our numerical results we compare them to corresponding results for non-fractal media to provide a context for their interpretation. In particular, we discuss the behavior of the coefficients in the limit of very large times. This case is relevant, since for practical applications it is commonly assumed that the ensemble and effective dispersion coefficients both converge at large times to the same constant value, called macrodispersion coefficients. Our results will therefore be important for the validity of this approach when applied to fractal media.

In the last part, we draw conclusions for the appropriate macro-scale transport model that should be used for the different scenarios discussed herein.

3.1 Explicit form of the dispersion coefficients in fractal media

In fractal media we have the following form for the ensemble and the effective dispersion coefficients

$$\begin{aligned} \mathscr {\widehat{D}}_{jj}^{\mathrm{ens}} (t)= & {} \int \limits _0^{\infty } d \varLambda \left\{ \overline{u} \, \widehat{L}_{1} \, \widehat{M}_j^- \left( \widehat{T}_u, 0 \right) \right\} \end{aligned}$$
(16)
$$\begin{aligned} \mathscr {\widehat{D}}_{jj}^{\mathrm{eff}} (t)= & {} \mathscr {\widehat{D}}_{jj}^{\mathrm{ens}} (t) - \int \limits _0^{\infty } d \varLambda \left\{ \overline{u} \, \widehat{L}_{1} \, \widehat{M}_j^+ \left( \widehat{T}_u, \widehat{T}_D \right) \right\} \end{aligned}$$
(17)

in which the arguments of the auxiliary functions are defined as follows

$$\begin{aligned} \widehat{T}_u= & {} \frac{t}{\widehat{\tau }_u} = \sqrt{2} \frac{t \, \overline{u} \, \lambda }{l}, \end{aligned}$$
(18)
$$\begin{aligned} \widehat{T}_D= & {} 2 \frac{t}{\widehat{\tau }_D} = 4 \frac{t \, \lambda ^2 \, D}{l^2}. \end{aligned}$$
(19)

The auxiliary functions \(\widehat{M}_j^{\pm } = f\left( \widehat{T}_u, \widehat{T}_D \right) \) themselves are given in isotropic fractal media in the longitudinal direction by

$$\begin{aligned} \widehat{M}_1^{\pm }= & {} \mp \, q_{uu} \sqrt{\frac{\pi }{2}} \frac{1}{\left( 1+2\widehat{T}_D \right) ^2} \bigg \{ \mathrm{erf} \left( g(\mp \widehat{T}_u) \right) \nonumber \\&+ \frac{1}{\sqrt{\pi }} \exp \left( -g^2(\mp \widehat{T}_u) \right) \nonumber \\&\times \left( \frac{1}{g(\mp \widehat{T}_u)} + \frac{4 \widehat{\psi }^2 w(\widehat{T}_u)}{g^2(\mp \widehat{T}_u)} - \frac{3}{2 g^3(\mp \widehat{T}_u)} \right) + \mathrm{erf} \left( g(\mp \widehat{T}_u) \right) \nonumber \\&\times \left( \frac{4 \widehat{\psi }^2 w(\widehat{T}_u)}{g(\mp \widehat{T}_u)} + \frac{8 \widehat{\psi }^4 w(\widehat{T}_u)}{g(\mp \widehat{T}_u)} - \frac{1}{g^2(\mp \widehat{T}_u)}\right. \nonumber \\&\left. - \frac{2 \widehat{\psi }^2 w (\widehat{T}_u)}{g^3(\mp \widehat{T}_u)} + \frac{3}{4 g^4(\mp \widehat{T}_u)} \right) \nonumber \\&- 8 \widehat{\psi }^4 \exp \left( \frac{1}{2 \widehat{\psi }^2} \right) \left[ \mathrm{erfc} \left( \frac{1}{\sqrt{2} \widehat{\psi }} \right) - \mathrm{erfc} \left( w(\widehat{T}_u) \right) \right] \nonumber \\&- \frac{4 \sqrt{8} \widehat{\psi }}{\sqrt{ \pi }} \left( \frac{1}{3} + \widehat{\psi }^2 \right) \bigg \} \end{aligned}$$
(20)

and in the transverse direction by

$$\begin{aligned} \widehat{M}_{2,3}^{\pm }= & {} \mp \, q_{uu} \sqrt{\frac{\pi }{8}} \frac{1}{(1+2\widehat{T}_D)^2} \bigg \{ \frac{1}{\sqrt{2 \pi }} \exp \left( -g^2(\mp \widehat{T}_u) \right) \nonumber \\&\left( - \frac{8 \widehat{\psi }^2 \, w(\widehat{T}_u)}{g^2(\mp \widehat{T}_u)} + \frac{3}{g^3(\mp \widehat{T}_u)} \right) \nonumber \\&+\, \mathrm{erf} \left( g(\mp \widehat{T}_u) \right) \left( - \frac{2 \widehat{\psi }^2 w(\widehat{T}_u)}{g(\mp \widehat{T}_u)} - \frac{8 \widehat{\psi }^4 w(\widehat{T}_u) }{g(\mp \widehat{T}_u)}\right. \nonumber \\&\left. + \frac{1}{2 g^2(\mp \widehat{T}_u)} - \frac{2 \widehat{\psi }^2 w(\widehat{T}_u)}{g^3(\mp \widehat{T}_u)} - \frac{3}{4 g^4(\mp \widehat{T}_u)} \right) \nonumber \\&- \left( 8 \widehat{\psi }^4 - 2 \widehat{\psi }^2 \right) \exp \left( \frac{1}{2 \widehat{\psi }^2} \right) \left[ \mathrm{erfc} \left( w(\widehat{T}_u) \right) \right. \nonumber \\&\left. -\, \mathrm{erfc} \left( \frac{1}{\sqrt{2} \widehat{\psi }} \right) \right] \nonumber \\&+ \frac{ \sqrt{8} \widehat{\psi } }{\sqrt{ \pi }} \left( \frac{1}{3} + 4 \widehat{\psi }^2 \right) \bigg \}. \end{aligned}$$
(21)

Here, we used the following definitions

$$\begin{aligned} q_{uu} = \frac{\sigma ^2_Y}{\gamma ^2} u^2, \quad \quad \widehat{\tau }_u = \frac{\widehat{L}}{\bar{u}}, \quad \quad \widehat{\tau }_D = \frac{\widehat{L}^2}{D}, \quad \quad \widehat{\varepsilon } = \frac{\widehat{\tau }_u}{\widehat{\tau }_D} \end{aligned}$$
(22)

in which \(q_{uu}\) is the scaled mean velocity squared, \(\widehat{\tau }_u\) is the advective time scale for the distance \(\widehat{L}\), \(\widehat{\tau }_D\) is the corresponding local-dispersive time scale and \(\widehat{\varepsilon }\) is the inverse Péclet number. We also used the following short notations

$$\begin{aligned} w(\widehat{T}_u)= & {} \frac{1}{\sqrt{2}} \frac{ \left( 1+2\, \widehat{T}_D\right) \frac{1}{\widehat{\varepsilon }} + \frac{t}{\widehat{\tau }_u}}{\sqrt{1+2\, \widehat{T}_D + \frac{2 \, t}{\widehat{\tau }_D}}}, \end{aligned}$$
(23)
$$\begin{aligned} g(\mp \widehat{T}_u)= & {} \frac{1}{\sqrt{2}} \frac{ \mp \frac{t}{\widehat{\tau }_u}}{\sqrt{1+2\, \widehat{T}_D + \frac{2 \, t}{\widehat{\tau }_D}}}, \end{aligned}$$
(24)
$$\begin{aligned} \widehat{\psi }= & {} \frac{\widehat{\varepsilon }}{\sqrt{1+2\, \widehat{T}_D}}. \end{aligned}$$
(25)

The explicit expressions for the longitudinal and transverse components of \(\mathscr {\widehat{D}}^{\mathrm{ens}} (t)\) can be derived by inserting Eqs. (20) and (21) into Eq. (16) and of \(\mathscr {\widehat{D}}^\mathrm{eff} (t)\) by inserting Eqs. (20) and (21) into Eq. (17), respectively. Due to the length of the resulting expressions we present them in the appendix only.

3.2 Longitudinal dispersion

For non-fractal media, the longitudinal and transverse component of the dispersion coefficient are known to behave quite differently. For fractal media, we saw the same behavior and therefore will present them separately here. Let us begin with the longitudinal component of the ensemble and effective dispersion coefficients. To investigate their behavior, we evaluated the corresponding expression [see Eqs. (48) and (52) in the Appendix] numerically.

Fig. 1
figure 1

The longitudinal dispersion coefficients \(\widehat{\mathscr {D}}^{[\cdot ]}_{11}\) in fractal media for \(0< \beta < 1 \)

Full size image
Fig. 2
figure 2

The longitudinal dispersion coefficients \(\widehat{\mathscr {D}}^{[\cdot ]}_{11}\) in fractal media for \(1< \beta < 2 \)

Full size image

Our results showed that both in the fully-fractal (see Fig. 1) as well as the semi-fractal (see Fig. 2) regimes the behavior of the longitudinal ensemble and effective dispersion coefficients differed substantially. The ensemble coefficient was larger than the effective dispersion values at all times, with the amount of the difference increasing with the degree of fractality, i.e., lower \(\beta \) (see Figs. 1, 2). However, in the semi-fractal regime \( 1< \beta < 2 \) the difference between the scales of the two quantities was smaller compared to the fully-fractal regime of \( 0< \beta < 1\) leaving the two quantities in the same order of magnitude.

Asymptotic behavior The long-term temporal behavior of the longitudinal coefficient can be derived by performing an asymptotic analysis of Eqs. (48) and (52) for the ensemble and effective part respectively. For the fully-fractal regime we got the following behavior

$$\begin{aligned} \mathscr {\widehat{D}}_{11}^{\mathrm{ens}, \infty } (t)= & {} c \cdot t^{1 - \beta } \end{aligned}$$
(26)
$$\begin{aligned} \mathscr {\widehat{D}}_{11}^{\mathrm{eff}, \infty } (t)= & {} c \cdot t^\frac{1-\beta }{2}, \end{aligned}$$
(27)

with c being a hitherto unknown constant. Equations (26) and (27) show that in the fully fractal regime the longitudinal ensemble and effective quantities are different at all times and represent two different properties. This result confirms our numerical investigations (Figs. 1, 2). This observation can be explained with the different temporal behaviors of the two quantities. Both are characterized by different power law relations given by Eqs. 26 and 27 (consistent with the findings of Fiori 2001). The different exponents depend on the degree of fractality. The different exponent for the ensemble quantity indicates the significant difference of the values. Furthermore, the power law relation shows that the dispersion values increase with a higher degree of fractality.

Note that the asymptotics given by Eqs. (26) and (27) are consistent with the results of Fiori (2001). In addition to this purely approximate behavior, we could, however, also estimate the constants in both expressions, therefore extending these previous results. To that end, we approximated Eqs. (48) and (52) for the marginal case of high Péclet numbers, i.e., small local dispersion coefficients. Our numerical investigations showed that these approximations were close to the case of finite Péclet numbers as well (results not shown). However, such power-law behavior is indicative of non-Fickian transport, i.e., the classical advection–dispersion equation with constant coefficients can not describe such behavior even in the long run. As a result, the macro-scale transport model has to be revised to account for this behavior. The resulting new model will be presented below, where we also discuss how the equivalent effective and ensemble longitudinal dispersion coefficients must be adapted and how their description can be derived from the asymptotic behavior of Eqs. (48) and (52).

3.3 Transverse dispersion

Now let us turn to the transverse component of the dispersion tensor. Compared to longitudinal component, the transverse component showed a wider range of relevant properties. In the following, we therefore describe our results pertaining to the ensemble and effective dispersion coefficient separately.

3.3.1 Transverse ensemble dispersion

The results for the evolution of the transverse components of the effective dispersion coefficient were obtained by numerically evaluating Eq. (51).

Influence of the degree of fractality To investigate the influence of the fractality on the transverse ensemble dispersion coefficient \(\widehat{\mathscr {D}}_{22}^\mathrm{ens}\), we evaluated Eq. (51) with varying values of \(\beta \) for both regimes, i.e. \(0< \beta < 1\) and \(1< \beta < 2\).

Fig. 3
figure 3

The ensemble transverse dispersion component \(\widehat{\mathscr {D}}_{22}^{\mathrm{ens}} \) in fractal media for \(0< \beta < 1\)

Full size image
Fig. 4
figure 4

The ensemble transverse dispersion component \(\widehat{\mathscr {D}}_{22}^{\mathrm{ens}} \) in fractal media for \(1< \beta < 2\)

Full size image

Contrary to the case of the longitudinal dispersion coefficient \(\widehat{\mathscr {D}}_{11}^{[\cdot ]}\), we saw an initial increase for \(\widehat{\mathscr {D}}_{22}^{\mathrm{ens}}\) followed by a monotonically decreasing behavior at later times in both regimes (see Figs. 3, 4). The values of \(\widehat{\mathscr {D}}_{22}^{\mathrm{ens}}\) differed dramatically between both regimes, being generally about three orders of magnitude larger in the fully-fractal than in the semi-fractal regime. This peak, which is typical for the transverse ensemble dispersion coefficient, was visible in both regimes but was more pronounced in the fully-fractal regime. Another noticeable difference to \(\widehat{\mathscr {D}}_{11}^{[\cdot ]}\) was the fact that in the large-time limit a single asymptotic value was reached. This value too, was different depending on the regime.

Influence of the micro-length\(\widehat{l}\)on\(\mathscr {\widehat{D}}_{22}^{\mathrm{ens}} (t) \) To investigate the influence of the length scale \(\widehat{l}\) on the transverse ensemble dispersion coefficient \(\widehat{\mathscr {D}}_{22}^{\mathrm{ens}}\) we evaluated Eq. (51) with varying values of \(\widehat{l}\).

Fig. 5
figure 5

The influence of the length scale \(\widehat{l}\) on the transverse ensemble coefficient \(\widehat{\mathscr {D}}_{22}^\mathrm{ens} \) for \(\beta = 0.5\)

Full size image

In our investigations, the influence of the length \(\widehat{l}\) on the transverse ensemble dispersion coefficient was very strong (see Fig. 5). For larger \(\widehat{l}\), the dispersion coefficient increased due to longer pathways that solute particles undergo and a higher spreading caused by a high variability in the flow field. Compared with the temporal behavior caused by a higher degree of fractality (see Fig. 3), the dispersion values for highly fractal media (\(\beta = 0.3\)) had almost the same magnitude as the ones observed for large lengths (\(\widehat{l} = 5\)) indicating long flow paths of solute particles. Furthermore, it could be seen that the length scale \(\widehat{l}\) had a more significant influence on the transverse ensemble coefficient. It affected not only the height of the peak but also its tailing.

Asymptotic behavior With increasing time, the transverse ensemble dispersion coefficient reaches an asymptotic limit which depends on the local dispersion

$$\begin{aligned} \lim _{t\rightarrow \infty } \mathscr {\widehat{D}}_{22}^{\mathrm{ens}} (t)= & {} \int \limits _0^{\infty } d \varLambda \bigg \{q_{uu} \, \overline{u} \, \frac{l}{\sqrt{2} \, \lambda } \sqrt{\frac{\pi }{2}} \, \bigg [\frac{2}{3 \sqrt{\pi }} \widehat{B} - 2 \widehat{B}^2\nonumber \\&+ \frac{16}{\sqrt{\pi }} \widehat{B}^3 -16 \widehat{B}^4 \nonumber \\&+ \exp \left( \frac{1}{4 \, \widehat{B}^2} \right) \mathrm{erfc} \left( \frac{1}{\sqrt{4 \, \widehat{B}^2}} \right) \nonumber \\&\left( 16 \, \widehat{B}^4 -2 \, \widehat{B}^2\right) \bigg ] \bigg \} \end{aligned}$$
(28)

with \(q_{uu}\) as given by Eq. (22) and

$$\begin{aligned} \widehat{B} = \frac{D \lambda }{l \bar{u}} \end{aligned}$$
(29)

being related to the timescale of the local dispersion. This behavior is similar to the transverse dispersion component in non-fractal media (Dentz 2000), but the time at which the asymptotic limit is reached is much larger in fractal media. It depends on \(\beta \) such that smaller values of \(\beta \) lead to larger values of the asymptotic limit as well as it being later reached. This analytical result therefore substantiates our numerical investigations presented above (see Fig. 3).

3.3.2 Transverse effective dispersion

The results for the evolution of the transverse components of the effective dispersion coefficient were obtained by numerically evaluating Eq. (55).

Influence of the degree of fractality To investigate the influence of the fractality on the transverse effective dispersion coefficient \(\widehat{\mathscr {D}}_{22}^\mathrm{eff}\), we evaluated Eq. (55) with varying values of \(\beta \) for both regimes, i.e. \(0< \beta < 1\) and \(1< \beta < 2\).

Fig. 6
figure 6

The effective transverse component \(\widehat{\mathscr {D}}_{22}^{\mathrm{eff}} \) in fully fractal media, i.e. for \( 0< \beta < 1 \)

Full size image
Fig. 7
figure 7

The effective transverse component \(\widehat{\mathscr {D}}_{22}^{\mathrm{eff}} \) in fully fractal media, i.e. for \( 1< \beta < 2 \)

Full size image

Our results showed that values of the effective transverse dispersion coefficient \(\widehat{\mathscr {D}}_{22}^{\mathrm{eff}}\) are generally larger in the fully-fractal regime compared to the semi-fractal regime (see Figs. 6, 7). In addition, the latter regime showed almost no influence of the fractality \(\beta \) on \(\widehat{\mathscr {D}}_{22}^{\mathrm{eff}}\).

Influence of the micro-length\(\widehat{l}\)on\(\mathscr {\widehat{D}}_{22}^{\mathrm{eff}} (t)\) To investigate the influence of the length \(\widehat{l}\) on the transverse effective dispersion coefficient \(\widehat{\mathscr {D}}_{22}^{\mathrm{eff}}\), we evaluated Eq. (55) with varying values of \(\widehat{l}\).

Fig. 8
figure 8

The influence of the length scale \(\widehat{l}\) on the transverse effective coefficient for \(\beta = 0.5\)

Full size image

Our results showed only a marginal effect of \(\widehat{l}\) on \(\mathscr {\widehat{D}}_{22}^{\mathrm{eff}}\) (see Fig. 8). This behavior can be explained with the fact that the correlation length (containing the length \(\widehat{l}\)) is not an important factor for the transverse effective dispersion coefficient, when compared to the local dispersion coefficient D (see Fig. 9 below).

Influence of the local dispersion coefficientDon\(\mathscr {\widehat{D}}_{22}^{\mathrm{eff}} (t) \) To investigate the influence of the local dispersion coefficient D on the transverse effective dispersion coefficient \(\widehat{\mathscr {D}}_{22}^{\mathrm{eff}}\), we evaluated Eq. (55) with varying values of D.

Fig. 9
figure 9

The influence of the local dispersion D on the transverse effective coefficient \(\mathscr {\widehat{D}}_{22}^\mathrm{eff}\) for \(\beta = 0.5\)

Full size image

Our results showed a significant impact of D on \(\widehat{\mathscr {D}}_{22}^{\mathrm{eff}}\) (see Fig. 9). The larger D, the larger was the transverse effective dispersion coefficient \(\widehat{\mathscr {D}}_{22}^{\mathrm{eff}}\). Compared with the effect of the length \(\widehat{l}\) (see Fig. 8), the local dispersion coefficient had a stronger impact on the transverse effective dispersion coefficient.

Asymptotic behavior At larger times, the transverse component of the effective dispersion tensor eventually reaches an asymptotic limit depending on the local dispersion coefficient

$$\begin{aligned} \lim _{t\rightarrow \infty } \mathscr {\widehat{D}}_{22}^{\mathrm{eff}} (t)= & {} \int \limits _0^{\infty } d \varLambda \bigg \{q_{uu} \overline{u} \frac{l}{\sqrt{2} \lambda } \sqrt{\frac{\pi }{2}} \, \bigg [\frac{2}{3 \, \sqrt{\pi }} \widehat{B}\nonumber \\&- 2 \widehat{B}^2 + \frac{16}{\sqrt{\pi }} \widehat{B}^3 -16 \widehat{B}^4 \nonumber \\&+ \exp \left( \frac{1}{4 \, \widehat{B}^2} \right) \mathrm{erfc} \left( \frac{1}{\sqrt{4 \, \widehat{B}^2}} \right) \nonumber \\&\left( 16 \, \widehat{B}^4 -2 \, \widehat{B}^2\right) \bigg ] \bigg \} \end{aligned}$$
(30)

with \(\widehat{B}\) given in Eq. (29) and related to the dispersive timescale. In the limiting case of no local dispersion (\(D = 0 \) and \(\widehat{\varepsilon } = 0\)), the asymptotic limit decreases to

$$\begin{aligned} \lim _{D\rightarrow 0} \lim _{t\rightarrow \infty } \mathscr {\widehat{D}}_{22}^{\mathrm{eff}} (t) = 0. \end{aligned}$$
(31)

Equation (30), describing the asymptotic value of \(\mathscr {\widehat{D}}_{22}^{\mathrm{eff}} (t)\) and Eq. (28), describing the asymptotic value of \(\mathscr {\widehat{D}}_{22}^{\mathrm{ens}} (t)\), are identical and therefore result in the same value. Hence we could conclude that

$$\begin{aligned} \mathscr {\widehat{D}}_{22}^{\infty } \equiv \lim _{t\rightarrow \infty } \mathscr {\widehat{D}}_{22}^{\mathrm{ens}} (t) = \lim _{t\rightarrow \infty } \mathscr {\widehat{D}}_{22}^{\mathrm{eff}} (t). \end{aligned}$$
(32)

This means that in the transverse direction an ergodic situation can be found in fractal media, too (compare Figs. 3, 6).

3.3.3 Comparing transverse ensemble and effective dispersion

So far we have seen that, contrary to the longitudinal dispersion coefficients, the transverse dispersion coefficients both reached an asymptotic value after a transient period. Moreover, Eq. (32) revealed that this asymptotic value is the same for the ensemble and effective quantity.

Fig. 10
figure 10

Asymptotic behavior of the transverse components \(\widehat{\mathscr {D}}_{22} \) in fully fractal media, i.e. for \( 0< \beta < 1 \)

Full size image
Fig. 11
figure 11

Asymptotic behavior of the transverse components \(\widehat{\mathscr {D}}_{22} \) in fully fractal media, i.e. for \( 1< \beta < 2 \)

Full size image

During the transient period however, strong differences could be seen between both quantities (see Figs. 10, 11). Regardless of the fractal regime, the transient period was always much longer (several orders of magnitude) for the effective quantity compared to the ensemble quantity. This difference was even more pronounced in the fully-fractal regime compared to the semi-fractal regime.

3.4 Comparing dispersion in fractal and non-fractal media

Following, e.g., Di Federico and Neuman (1998a), we compared our results for fractal media to results for non-fractal media which have been previously reported (Dentz et al. 2000a, 2002). This will help to provide a context for our results by highlighting differences and similarities between fractal and non-fractal media.

3.4.1 Longitudinal dispersion

Longitudinal dispersion can be conceptualized such that particles entering a medium at the same time may leave at different times due to the different flow paths they follow. The nature of these flow paths is therefore crucial for the observed longitudinal dispersion Fig. 12.

Fig. 12
figure 12

The longitudinal components of the ensemble and effective dispersion coefficients in non-fractal media according to Dentz et al. (2000a, 2002)

Full size image

In heterogeneous media, particles may experience different travel speeds, and therefore different mean square displacement, depending on the flow path they randomly follow. In non-fractal media, this difference vanishes with time since all heterogeneities will have been sampled equally by the particles after a given time. After a certain point in time, the mean travel speed of all particles therefore becomes identical regardless of their corresponding flow path. As a result, the longitudinal dispersion grows during an initial transient period until it reaches an asymptotic limit (Dentz et al. 2000a).

But this is not necessarily the case in fractal media, where a particle may perpetually experience new and stronger heterogeneities along its path. The occurrence of such new extreme values may become rarer over time but does not necessarily become zero. The longitudinal dispersion coefficient may therefore grow indefinitely as well. This situation was indeed observed in the fully-fractal regime (see Fig. 2), where no asymptotic limit was reached and the ensemble and effective quantities were characterized by different power laws at all times. In the semi-fractal regime, the two quantities eventually reached an asymptotic limit, but it was different for the two quantities (Fig. 2). These observations imply that ensemble and effective dispersion indeed describe two different transport properties. The timescale, at which the asymptotic limit is reached, is much shorter in non-fractal media because the flow paths of the solute particles are smoother and thus not as long as in fractal media. It is worthwhile to note that this stark difference between the fractal and non-fractal behavior was not observed by a similar comparison of Di Federico and Neuman (1998a). The reason for this discrepancy is simply the aforementioned notion that particles in fractal media may indefinitely experience new heterogeneities. In a truncated fractal medium, as assumed by Di Federico and Neuman (1998a), this is no longer the case once the mean travel length of the particle has exceeded the upper truncation length. After this point, the particles would have sampled all heterogeneity and the transport behavior would become essentially identical to that in a classic non-fractal medium.

3.4.2 Transverse dispersion

Transverse effective dispersion can be conceptualized such that particles entering a medium at the same point may leave at different points, due to local spreading (molecular diffusion at the pore scale and local dispersion at the continuum scale) causing the particles to randomly jump to different flow paths. On the other hand, transverse ensemble dispersion can be conceptualized as spreading caused by the different realizations of the porous medium (Fig. 13).

Fig. 13
figure 13

The transverse components of the ensemble and effective dispersion coefficients in non-fractal media according to Dentz et al. (2000a, 2002)

Full size image

In non-fractal heterogeneous media, the transverse dispersion coefficient grows initially, which is in analogy to the observed longitudinal dispersion. Eventually, the value of the transverse dispersion coefficient reaches a peak after which it decreases again, finally saturating at an asymptotic limit. Compared to the longitudinal case, this difference in behavior was attributed by Dentz et al. (2002) to the incompressibility of the fluid, whereas Sposito (2001) connects the same behavior to the topological properties of the flow field, namely having zero helicity density everywhere. Dispersion is impacted by the effective length of flow paths and the velocity experienced along them. In both the longitudinal and transverse direction, the heterogeneities cause the dispersion coefficient to grow but in the latter, the properties of the flow field prevent the streamlines from diverging too far and leading to a comparably poor mixing overall.

Since we saw qualitatively the same behavior in fractal media, we conclude that the same properties of the flow field are dominating in fractal media, too. The main difference between the two types of media was the scale of transverse dispersion values and time. In addition, the asymptotic value in fractal media was higher for smaller values of \(\beta \), i.e. a higher fractality, and the asymptotic limit was reached at a later point of time (Figs. 10, 11 ).

3.5 Asymptotic transport model for fractal media

At this point, we have seen that in fractal media, too, an asymptotic limit is reached in the semi-fractal regime, i.e., \(1< \beta < 2\). The initial time period, during which the behavior is non-Fickian is, however, much longer compared to the non-fractal case. Using the asymptotic values provided above with a classical ADE [see Eq. (1)] for simulation purposes is therefore possible but less reliable compared to the non-fractal case.

In addition, we saw the non-Fickian behavior of longitudinal dispersion indefinitely prevailing in the fully-fractal regime, i.e., \(0< \beta < 1\). Modeling the transport behavior with a classical ADE is therefore always bound to produce errors, a notion that has already been acknowledged by earlier studies (Niemi et al. 2000; Liu et al. 2007; Reeves et al. 2008; Roubinet et al. 2012). To account for this discrepancy, we propose the following asymptotic model

$$\begin{aligned}&\frac{\partial }{\partial t} c^{[\cdot ]}( {\mathbf {x}}, t)\nonumber \\&\quad = \nabla \cdot \left( \begin{pmatrix} \mathscr {D}^{[\cdot ]}_{11} &{} 0 &{} 0 \\ 0 &{} \mathscr {\widehat{D}}^{\infty }_{22} &{} 0 \\ 0 &{} 0 &{} \mathscr {\widehat{D}}^{\infty }_{22} \end{pmatrix} \begin{pmatrix} \frac{\partial ^{\alpha -1}}{\partial |x|^{\alpha -1}} \\ \frac{\partial }{\partial y} \\ \frac{\partial }{\partial z} \end{pmatrix} - {\mathbf {u (x)}} \right) c^{[\cdot ]}({\mathbf {x}}, t). \end{aligned}$$
(33)

In the fully-fractal regime the resulting asymptotic model is therefore still a classical ADE in the transverse directions but has been revised to a fADE in the longitudinal direction (Chen and Hsu 2007; Zhang et al. 2008; Hsu and Chen 2010). The choice of the coefficient \(\alpha \) in Eq. (33) is following the established notation such that a value of \(\alpha = 2\) would recover a regular Gaussian transport regime (Schumer et al. 2009).

Compared to a classical ADE, a major difference is the relationship between the coefficients described so far and the corresponding coefficients \(\mathscr {D}^{[\cdot ]}_{11}\) of the governing differential equation, as given in Eq. (33). In the classical case, both are identical, whereas in the case of fractal transport behavior, the relationship is not as simple. This can be demonstrated by the corresponding dimensions, with fractional dispersion coefficients having a fractal dimension of \(L^{\alpha } T^{-1}\), which has no correspondence in the coefficients describing the behavior of the plume. To derive the former, we connect the coefficient of the fADE \(\mathscr {D}^{[\cdot ]}_{11}\) and the coefficient describing the plume behavior \(\mathscr {\widehat{D}}^{[\cdot ]}_{11}\) through the spread of the plume

$$\begin{aligned} \sigma _{11}^2 = (\mathscr {D}^{[\cdot ]}_{11}\ t)^{\frac{2}{\alpha }} = \mathscr {\widehat{D}}^{[\cdot ]}_{11}\ t. \end{aligned}$$

This expression is motivated by matching the dimensions of each dispersion coefficient to \(\sigma _{11}^2\). Consequently, we can derive an expression for \(\mathscr {D}^{[\cdot ]}_{11}\) by rearranging this expression, i.e.

$$\begin{aligned} \mathscr {D}^{[\cdot ]}_{11}= & {} (\mathscr {\widehat{D}}^{[\cdot ]}_{11}\ t)^{\frac{\alpha }{2}} t^{-1} \\= & {} ({\mathscr {\widehat{D}}^{[\cdot ]}_{11}})^{\frac{\alpha }{2}} t^{\frac{\alpha }{2}-1} \end{aligned}$$

Inserting the asymptotic expression for the dispersion coefficients \(\mathscr {\widehat{D}}^{[\cdot ]}_{11}\), would therefore yield the expression for the respective \(\mathscr {D}^{[\cdot ]}_{11}\) in the asymptotic transport model.

Ensemble dispersion Using these statements, the ensemble spreading of the solute cloud is evolving for large t asymptotically with

$$\begin{aligned} \mathscr {\widehat{D}}^{\mathrm{{ens}}}_{11} = \frac{\sigma ^2_Y}{1-\beta } \left( \frac{ut}{L} \right) ^{1-\beta } u L. \end{aligned}$$

Please note that in this expression, the \(\sigma ^2_Y\) is referring to the variance of the log-conductivity field Y and not to the spread of the solute cloud. Inserting this expression into above relationship for the asymptotic behavior of \(\mathscr {D}^{[\cdot ]}_{11}\), we get

$$\begin{aligned} \mathscr {D}^{\mathrm{{ens}}}_{11} = \left( \frac{\sigma ^2_Y}{1-\beta } \left( \frac{ut}{L} \right) ^{1-\beta } u L \right) ^{\frac{\alpha ^{\mathrm{{ens}}}}{2}} t^{\frac{\alpha ^{\mathrm{{ens}}}}{2}-1} \quad for \quad 0< \beta < 1 \end{aligned}$$
(34)

with

$$\begin{aligned} \alpha ^{\mathrm{{ens}}} = \frac{2}{2 - \beta }. \end{aligned}$$

This last expression is derived by comparing the power laws of Eqs. (34) and (26).

Effective dispersion Using again the statements above, the expected spreading of a single solute cloud is evolving for large t asymptotically with

$$\begin{aligned} \mathscr {\widehat{D}}^{\mathrm{{eff}}}_{11} = \sigma ^2_Y u L \left( \frac{D}{L^2} t \right) ^{\frac{1 - \beta }{2}} f(\beta ) \end{aligned}$$

with

$$\begin{aligned} f(\beta ) = \frac{\sqrt{2}}{\varGamma (\beta /2)} \int d\lambda \lambda ^{\beta - 1} \left( 1 - (1 + 8 \lambda ^2)^{-1/2} \right) . \end{aligned}$$

This results in

$$\begin{aligned} \mathscr {D}^{\mathrm{{eff}}}_{11} = \left( \sigma ^2_Y u L \left( \frac{D}{L} t \right) ^{\frac{1 - \beta }{2}} f(\beta ) \right) ^{\frac{\alpha ^{\mathrm{{eff}}}}{2}} t^{\frac{\alpha ^{\mathrm{{eff}}}}{2}-1} \quad for \quad 0< \beta < 1 \end{aligned}$$
(35)

with

$$\begin{aligned} \alpha ^{\mathrm{{eff}}} = \frac{4}{3 - \beta } \end{aligned}$$

This last expression is again derived by comparing the power laws of Eqs. (35) and (27).

Equations (34) and (35) now provide the expressions that can be used to specify the fractional dispersion coefficient. Inserting either expression into Eq. (33) will complete the asymptotic transport model for \(c^{\mathrm{{ens}}}( {\mathbf {x}}, t)\) and \(c^{\mathrm{{eff}}}( {\mathbf {x}}, t)\) respectively.

4 Critical assessment of the transport model

Having now derived the effective transport model for three-dimensional, isotropic media, we would like to finally assess its applicability and limitations vis-à-vis the assumptions being made here as well as other approaches from the literature.

Let us begin with the asymptotic transport model proposed here. Since this is an fADE in the longitudinal direction and a classic ADE in the transverse direction, both derived with a first-order analysis, a number of limitations have to be expected. First, like e.g., MRMT or CTRW, a fADE assumes a statistically stationary velocity field to underlie the transport behavior (Neuman and Tartakovsky 2009). This is problematic if, say, measurements of the groundwater velocity are available, since this incorporating these measurements makes the velocity field statistically instationary. Second, while the solute cloud does not need to be Gaussian, fADE can only generate symmetric solute clouds. However, depending on the heterogeneity of the conductivity field, asymmetrical solute clouds are possible and MRMT and CTRW models can generally reproduce asymmetric transport behavior. The validity of our transport model is therefore dependent on the validity of the used conductivity model. Furthermore, fADE as well as CTRW are often praised for their lack of a support scale, making them completely scale independent. Yet by explicitly including a local scale in our conductivity model, we violate this feature in our study. However, following, e.g., Neuman and Tartakovsky (2009), we consider this to be a possible benefit instead of liability since it allows us to connect the parametrization of our model to observed data of the conductivity field. Finally, fADE are sometimes criticized for failing to provide a link between the fractional parameters and measurable medium properties (Neuman and Tartakovsky 2009). In our study, we close this problematic gap by showing how \(\alpha \), the parameter of the fractional derivative, depends on \(\beta \), the parameter of the hydraulic conductivity field.

Now let us discuss the main theoretical assumption of this study, i.e. the fractal nature of the heterogeneous medium with an ever evolving scale of heterogeneities. Such a mathematical model has, however, limits for any real-world physical system. In particular, physical fractals are always constrained by upper and lower limits. Here, the lower limit is the small-scale resolution, i.e., the support scale, of the system as already discussed above. However, any real geological medium is also constrained by an upper length scale, i.e., its finite size. This means that the effective description, described herein, may be less reliable in situations where the solute cloud reaches the confines of the physical domain being investigated. While such finite-size effects are present in many effective macro-scale descriptions, its potential error-inducing effects must be acknowledged and proper care by the practitioner is advised.

Finally, let us turn to one of the main findings of this paper, namely that in fractal media, too, the transverse dispersion coefficient does not grow indefinitely. This was tentatively connected to either the incompressibility of the medium or topological constraints as discussed by Sposito (1994, 2001). However, recent studies have provided evidence that these topological constraints may not hold for all types of heterogeneous media (Cirpka et al. 2015; DiDato et al. 2016). While it is not clear whether such media do indeed increase transverse dispersion, it is prudent to say that the upper limit on transverse dispersion, reported herein, may not apply to all types of heterogeneous media.

5 Conclusions

For this study, we followed Neuman (1995) who proposed the existence of log-hydraulic conductivity fields with fractal properties. This fractal behavior causes the longitudinal ensemble and effective dispersion coefficients to increase with travel distance as the solute samples more and more heterogeneities. The temporal behavior of the longitudinal dispersion coefficient in such fractal media stands in contrast to classic conceptualizations of heterogeneous media, where the different scales are clearly separated. In classic, non-fractal media the ensemble and effective dispersion coefficients reach the same asymptotic value at large times. In contrast, in fractal media these quantities do not reach an asymptotic limit since the longitudinal dispersion coefficient is dominated by the infinite correlation length \(\widehat{L}\). The ensemble and effective quantities are different at all times because they are characterized by different power laws. This implies that they describe two different transport properties and for the longitudinal component we find a non-ergodic temporal behavior.

The behavior of the transverse ensemble and effective dispersion coefficients is comparable to that in classic, non-fractal media as shown, e.g., by Dagan (1988), Kitanidis (1988) and Musuuza et al. (2011) to first increase and then reduce with travel distance. However, a remarkable difference between the temporal behavior in the transverse dispersion coefficients in both media with clearly separated scales and fractal media with evolving scales is the much larger timescale to approach the asymptotic behavior in the latter. This large timescale is caused by the fractal and therefore much longer pathways of the solute particles. We showed how the timescale increases with the degree of fractality.

Having presented the behavior of the longitudinal and transverse dispersion coefficients, we analyze their asymptotic behavior. Using this analysis, we are able to derive an effective transport model for larger times. This model is a classic ADE in the transverse direction and a fADE in the longitudinal direction, accounting for the effects of the fractal medium. We present the coefficients of this model and show how they are connected to the fractal properties of the conductivity model.