The effects of heterogeneity on solute transport in porous media: anomalous dispersion

Mixing of solute in sub-surface flows, for example in the leakage of contaminants into groundwater, is quantified by a dispersion coefficient that depends on the dispersivity of the medium and the transport velocity. Many previous models of solute transport data have assumed Fickian behaviour with constant dispersivity, but found that the inferred dispersivity increases with the distance from the point of injection. This approach assumes that the dispersion on either side of the advective front is symmetric and that the concentration there is close to half of the injected value. However, field data show consistent asymmetry about the advective front. Here, motivated by experimental data and a fractal interpretation of porous media, we consider a simplified heterogeneous medium described by a dispersivity with a power-law dependence on the downstream distance from the source and explore the nature of the asymmetry obtained in the solute transport. In a heterogeneous medium of this type, we show that asymmetry in solute transport gradually increases with the increase in heterogeneity or fractal dimension, and the concentration at the advective front becomes increasingly different from 50% of that at the inlet. In particular, for a sufficiently heterogeneous medium, the concentration profiles at late-time approach a non-trivial steady solution and, as a result, the concentration at a downstream location will never reach that at the inlet. By fitting the Fickian solution to our results, we are able to connect the parameters from our model to those found from experimental data, providing a more physically grounded approach for interpreting them.


Introduction
Dispersion in porous media is the process whereby the tortuosity of flowpaths enhances the molecular mixing of constituents.It is therefore fundamental to the mixing and distribution of tracers -for example in the dispersal of contaminants in groundwater [1], the intrusion of seawater into groundwater [2], or the dispersion of geochemical tracers [3] -and has been studied extensively [4].Theoretical models of dispersion are generally based on an effective advection-diffusion equation, where C is the concentration, t is time, and x is the downstream distance from the source.U is the average flow velocity through the porous medium, related to the Darcy flux q and porosity φ via q = φU , while D is the effective diffusivity that incorporates both molecular diffusion and dispersivity.This approach follows the classic work by Taylor [5] on flow and dispersion in a pipe (conceptually a simplified porous medium), where dispersion is a consequence of the longitudinal stretching of the tracer through transverse shear.This results in greatly enhanced transverse diffusion rates which may be interpreted, at the pore scale, as enhanced longitudinal dispersion.In Taylor's work, transport of an initial step-function of a non-reactive tracer results in a concentration profile that may be described by self-similar functions that are symmetrical about the mean concentration, which travels at the mean fluid velocity.Taylor found that the effective dispersion coefficient in the direction of flow, may be related to the square of the fluid velocity at the axis u 0 , the pipe radius a, and the molecular diffusion coefficient D m .Since Taylor's original work, similar models of velocity-dependent dispersion have been applied in a number of contexts.In general, these studies find that compilations of laboratory measurements and field observations of the dispersion of tracers at a fixed distance from the source in porous media may be fit to the advection-dispersion equation assuming a constant effective dispersivity expressed as the sum of the molecular diffusion coefficient D m with the product of a longitudinal dispersivity α and the average flow velocity U [1,6].In many settings, it is assumed that the Péclet number Pe = αU /D m 1 such that the effects of molecular diffusion are negligible and so D ≈ αU . (4) However, experimental measurements show that the apparent 'constant' α increases with the distance travelled by the tracer.Neuman [7] suggests from a compilation of laboratory and field data that α scales with advective distance, x = U t, as where α 0 is a reference dispersivity and has dimensions of [L] 1−γ , and found that the observations were best explained by γ ≈ 1.5.After analysing the data from other field experiments, Gelhar et al. [6] disputed this, suggesting that there is no single relationship such as equation ( 5) for all porous media, and that the degree of heterogeneity in an aquifer is likely to control the magnitude and scale dependence of its dispersivity.
There have been numerous attempts to create theoretical models of porous flow and predict dispersivities as a function of the heterogeneity of the permeability and pore structure of the porous medium.Pickens and Grisak [8], following Mercado [9], showed that macro-dispersivity of an aquifer with horizontally layered permeability in the direction of flow and no exchange between layers would evolve as where K is the mean hydraulic conductivity and σ K is the standard deviation of the hydraulic conductivity.Thus, they infer a value of γ = 1.However, more complex models describing aquifer structure either by a fractal hierarchy of permeabilities [10] or as a stochastic distribution of permeabilities [11] predict a range of behaviours.In some models,the longitudinal dispersivity approaches a constant value at long travel times and at large travel distances, while at early times and at short distances, it exhibits 'non-Fickian' behaviour described by equation ( 5) with γ ≈ 1 [11,12].Pickens and Grisak [8], using multiple depth samplers in field experiments in a sand aquifer, observed that dispersion in individual homogeneous horizons was not scale-dependent but in the whole aquifer longitudinal dispersion scaled with distance (γ ≈ 1).The central assumption for most of these models is that dispersion occurs symmetrically on both sides of the mean advective distance, x, and that the relative concentration at x = x is half of the injected value at x = 0.However, a few researchers, e.g.Day [13], recognised that there is persistent asymmetry in the time-concentration curves of field data.The asymmetry in the data implies that the dispersion coefficient is not constant on either side of x.Furthermore, Gelhar et al. [1] cast doubt on the applicability of a Taylor-type dispersion model at the field scale and Pickens and Grisak [8] agree that there is no common consensus on the dependence of dispersivity on the advective distance and that other mathematical forms may be more applicable than what they present.
In recent decades, numerical techniques for modelling anomalous dispersion have gained tremendous popularity [14,15].Some of the popular methods use random realisations [16] or pore network models [17] to predict the behaviour of solute transport stochastically in a representative elementary volume.Given that aquifers are highly heterogeneous and that their boundaries are not well defined, techniques involving space-time-nonlocal representations are widely used [14].Some of the popular nonlocal models are a continuous time random walk model [18], fractional advection-dispersion models [19][20][21], and a Lagrangian model [22], among others.Some models include multirate mass transfer [23] and delayed diffusion techniques [24].All the above models have been successful to varying degrees in modelling solute transport observed in field settings.
These previous studies have sought to develop the physical picture underpinning the large scale effects of dispersion.Most of them, however, are indeterministic in nature.In contrast, in the present work, we present simple analytical results that we obtain using a power-law spatial dispersivity which is motivated by the concept of fractal hierarchy of heterogeneity, introduced by Wheatcraft and Tyler [10].The fractal hierarchy states that the degree of irregularity or the increase in heterogeneity in aquifers is self-similar and independent of scale.Building on this concept, here we assume that the self-similarity in the fractal heterogeneity applies to the dispersivity as well such that the increase in the value of dispersivity with the downstream distance is also self-similar.Therefore, we approximate the dispersivity by where α 0 is the reference dispersivity at x = 0 and has a dimensional unit of [L] 1−β , while β encodes the fractal dimension and the magnitude and correlation length of heterogeneity.In equation (7), dispersivity is assumed to depend on the actual downstream distance, x, away from the source rather than the advective distance x, in contrast to the previous consideration in equation (5).We thereby explore the asymmetry in the solutions to the advection-diffusion equation.
There have been a large number of studies of the advection-dispersion equation with a similar dispersivity to (7).Several have considered linear power-law dispersivity (i.e.β = 1) in one or more spatial dimension [25][26][27][28][29][30].Some have extended this to more complicated functional forms that are linear near the injection point and tend to a constant in the far field [31][32][33] or have considered dispersivities that scale with particular powers of x and t in order to derive self-similar solutions [34,35].Others have looked at spatial power-law dispersivity as in (7) [36,37], but have not described the full range of behaviours as in the present work.
In equation (7), the value of β = 0 indicates that the medium is homogeneous, whereas β > 0 implies a heterogeneous medium, with larger values of β corresponding to larger heterogeneity in the medium.It should be noted that previous interpretations of the longitudinal dispersion from single laboratory and field experiments are based on assuming a constant macro-scale dispersivity at a fixed travel distance or time, i.e. implicitly β = 0.In this work, we investigate the effects of heterogeneity and fractal hierarchy, quantified in the form of β, on solute transport in a simplified medium with dispersivity described by (7).The solutions from the reduced model presented here characterise the different behaviours that occur for a dispersivity of this type and could be used as a computationally efficient way to assess the degree of heterogeneity present in a real-world system.
The remainder of the paper is structured as follows.We first set out the mathematical model of the problem being considered in Sect.2, for both one-dimensional and radial flow.In Sect.3, we present the numerical solution of this problem, as well as the exact and approximate analytical solutions that are then discussed in more detail in Sect. 4. In Sect.5, we consider the breakthrough curves produced by this system, and the implications of interpreting it through a Fickian lens.We then discuss the broader implications of our results in Sect.6 and present our conclusions in Sect.7.

One-dimensional flow
We consider the one-dimensional flow of fluid with initial solute concentration C 0 into a porous medium saturated with ambient fluid with solute concentration C a .The average flow velocity of the injected fluid, U , is taken to be constant.The injected fluid is assumed to be completely miscible with the ambient fluid and therefore the conservation C of solute may be modelled by an effective advection-dispersion equation This is linear in C, so we can instead work with the reduced concentration The initial condition and boundary conditions are then We are interested here in the case when the Péclet number is large, so that molecular diffusivity can be neglected outside of a narrow region near the injection point, and the diffusivity due to dispersion scales as a power of the downstream distance x: Here, α 0 is a reference dispersivity and has dimensions of L 1−β .We will focus here on the behaviour of the solution for β ≥ 0, with β = 0 corresponding to the case of Fickian Diffusion.The full governing system for the concentration C(x, t) can then be written as We would expect the three terms in (12) to respectively scale like For the special case of β = 1, we can balance all three of these terms and look for a self-similar solution in terms of the similarity variable η = x/(U t).For β = 1, ( 16) instead suggests the rescaled variables This allows us to write the governing equations in the simpler form The expected scaling behaviours of the advective and dispersive terms now suggest For 0 ≤ β < 1, we thus expect the solution to be predominantly governed at early times by dispersion, when the main part of the solution lies in the region z < 1, and at late times by advection, when the solution has spread sufficiently far into the region z > 1. Correspondingly, for β > 1, we would expect the reverse: advection dominating at early times, and dispersion at late times.For the intermediate case of β = 1, as seen above, advection and dispersion remain in balance at all times.

Radial flow
Above we have set up the mathematical model for power-law dispersion in flow along one dimension.Equally important is the case of axisymmetric flow; thankfully, there is a convenient trick shown by Philip [36] that allows us to reduce the axisymmetric system to the one discussed above.
In the axisymmetric case, the average flow velocity now varies spatially as Ũ /r , due to conservation of mass, where we take Ũ to be constant.If we express the dispersion coefficient as D = α0 Ũr δ , then the governing advection-dispersion equation is now Following [36], we introduce y = (1/2)r 2 , so that (1/r )∂ r = ∂ y .In terms of y, (23) becomes Relabelling the parameters as we are left with (12).The parameter range −2 ≤ δ ≤ 0 corresponds to 0 ≤ β ≤ 1, and δ > 0 to β > 1, and so the results in the linear geometry can straightforwardly be applied to the case of axisymmetric flows.

Numerical results
The governing system (18)-( 21) was solved numerically via Finite Differences, using the Il'in upwinding discretisation scheme [38,39] to accurately capture both advection and dispersion.An adaptive, predictor-corrector time-stepping routine was used to keep the relative error due to a finite step size beneath a set relative threshold (typically 10 −7 ).The far-field boundary condition (20) was replaced by a zero gradient condition at the end of the numerical domain, z = z ∞ .An adaptive, non-uniform spatial grid was employed to minimise the effect of this artificial far-field boundary condition, with z ∞ doubled while the number of grid points was kept the same whenever C(z ∞ , T ) rose above a threshold value (typically 10 −7 ), while also maintaining resolution near the inlet and up to the advective front.The agreement shown in Fig. 1 between the numerically determined profiles and their corresponding exact or approximate analytic solutions reassures us that our solver has accurately solved this system.
Figure 1 highlights the full range of behaviour of this system for β ≥ 0. Plotted are the numerical results for β = 0, 0.5, 1, and 1.5 from early to late times, first in the advective frame in the left-hand column, and then against log(z) in the right-hand column.The analytical solutions discussed in Sect. 4 are also shown for comparison.
Figure 1(a) shows the Fickian case of β = 0, with the magenta dashed curves showing the known exact solution discussed in Sect.4.1; this is plotted at early and late times, but matches the numerical solution at all intermediate times as well.The two panels of Fig. 1(a) show how the profiles transition from a solution that is initially dominated by dispersion, due to the strong gradients in the initial step-function profile at T = 0, to one that is predominantly governed by advection at late times, with a dispersive region around the advective front that diminishes in time in the advective frame z/T as this region only grows like T 1/2 .
Figure 1(b) shows the β = 0.5 case, which is representative of 0 < β < 1.This follows a similar pattern to the Fickian case -governed by dispersion at early times and advection at late times.However, due to the form of the dispersivity, D ∼ z β , in comparison with the Fickian case dispersion here is initially weaker, in the region z < 1 and so at early times, and stronger at late times, when the bulk of the profiles lies in the region z > 1.The magenta dashed curves here correspond to the early-time and late-time approximate solutions discussed in Sect.4.2.
Figure 1(c) shows the dividing case of β = 1.As discussed in Sect.2, this is the case in which the advective and dispersive terms can be in balance at all times, producing a self-similar solution -this is clear from the profiles plotted in the advective frame, in terms of η = z/T , where all of the profiles collapse onto a single curve.The rescaling of variables to z and T as defined in ( 17) is not possible unless α 0 = 1, and so this is the particular case shown in Fig. 1(c) so that it can be compared with the other values of β considered.The magenta dashed curve here shows the self-similar profile for α 0 = 1, which is discussed for general α 0 in Sect.4.3.
Finally, Fig. 1(d) shows the solution for β = 1.5, which is representative of β > 1.In contrast to 0 ≤ β < 1, advection now dominates at early times.At late times, however, there is a non-trivial, physically relevant steady solution that the profiles accumulate towards, as can be seen in the right-hand panel of Fig. 1(d).As a result, the Fig. 1 Concentration profiles for β = 0, 0.5, 1, 1.5 at a logarithmically spaced range of times from T = 10 −3 to T = 10 3 .In the left column, the profiles are plotted in the advective frame, against z/T = x/(Ut), while in the right column, the same profiles are plotted against log z.The magenta dashed curves are the early-time and late-time analytical results discussed in Sect. 4 concentration at any fixed location z > 0 will never reach the value at the inlet.The magenta dashed curves here show the early-time and late-time approximate solutions discussed in Sect.4.4.

ˇ= 0: Fickian Dispersion
For the well-known case of Fickian Dispersion, corresponding to β = 0, there is an exact solution [40] that can be derived for example via a Laplace Transform, or the convenient method discussed in [36].This exact solution is or in terms of x and t, where erfc is the complementary error function [41]: At late times, the second term in (27) can be neglected away from the inlet at x = 0.The size of the dispersive region around the advective front then grows as √ α 0 U t.

Early-time behaviour
When β < 1, the ratio of the advective and dispersive scalings, (22), suggests that dispersion will dominate at early times.The balance between the time derivative and the dispersive term in (18) then points towards changing variable from z to where the prefactor here is chosen for convenience later in the analysis.With this substitution, we can rewrite (18) as where we have neglected the term corresponding to ∂C/∂z on the left-hand side of (18).This can be rewritten as which, on applying the inlet and far-field boundary conditions ( 19) & ( 20), has the self-similar solution where U is the upper incomplete gamma function [41]: This solution is shown for a range of β values in Fig. 2(a).Based on the term neglected in order to derive this approximation, the next-order correction to the solution is

Late-time behaviour
At late times, advection now dominates.If we switch to a frame moving with the advective front, i.e. from z to ζ = z − T , then (18) becomes where in the last step, we have assumed that dispersion is the strongest in the region around the advective front, where |z/T − 1| 1.This then suggests the similarity variable where the prefactor has been chosen for convenience below.The governing equation can then be written as Applying the inlet and far-field boundary conditions, ( 19) & (20), this points us towards the self-similar solution The normalisation constant N is chosen here so that C(0, T ) = 1: This solution is shown in Fig. 2(b).In comparison with the Fickian solution (27), the size of the dispersive region around the advective front here grows as α (38) agrees with the late-time behaviour of the Fickian solution (26).

ˇ= 1: Self-similar behaviour
As seen in Sect.2, the advection and dispersion terms in (12) are in balance at all times for β = 1.This allows us to seek a self-similar solution in terms of η = x/U t, which in fact turns out to be an exact solution for our initial and boundary conditions ( 13)-( 15) and has been studied previously [30,34].For C(x, t) = C(η), (12) becomes Integrating this twice and applying the inlet and far-field boundary conditions, ( 13) & ( 14), leads us to the solution When α 0 = 1, this simplifies to C(x, t) = e −x/Ut .Figure 3 shows the behaviour of this solution for different values of α 0 .Unlike for β = 1, where the dispersion coefficient α 0 could be scaled out of the governing equation via z and T , here it plays a direct role in determining the character of the solution.

Early-time behaviour
At early times for β > 1, the solution is now controlled primarily by advection and we can reuse the approximate solution derived in Sect.4.2.2: where

Late-time behaviour
At late times, advection no longer dominates, but we are not able to apply the dispersion-dominated approximation derived in Sect.4.2.1 -U (χ , ν) diverges as ν → 0 when χ = (1 − β)/(2 − β) < 0. Instead, there is a physically relevant steady solution C s (z) which C(z, T ) approaches at late time for β > 1, as seen in Fig. 1(d).This steady solution is given by For 0 ≤ β ≤ 1, the only physically relevant steady solution is the trivial one, C s (z) = 1.For β > 1, we are able to apply both the inlet and far-field boundary conditions, ( 19) & ( 20), to get The behaviour of C s (z) for different values of β is shown in Fig. 4.This steady solution has an important implication for a system with β > 1.Whereas for 0 ≤ β ≤ 1 the concentration at any fixed location z will tend to 1 as T → ∞, for β > 1 the concentration will instead tend to C s (z) < 1 -the concentration downstream will never reach that at the inlet.We can approximate how the concentration profile C(z, T ) approaches this steady solution C s (z).If we write the concentration as C(z, T ) = C s (z)F(z, T ) then, after some manipulation, the governing equation ( 18) becomes When z is sufficiently large, z −(β−1) 1 and we can use the approximation coth u ∼ u −1 + (1/3)u + O(u 3 ) for u 1.For 1 < β < 2, a scaling balance again suggests that we should introduce the dispersive variable ν defined in (29) and expand F as F(z, T ) = F 0 (ν) + . . . .At leading order, we get In a similar manner to (32), we arrive at The second term in the expansion of coth u above suggests that the next-order correction to the solution is O(T −2(β−1)/(2−β) ).For β ≥ 2, a similar, but more involved, process can be followed to determine the late-time behaviour.

Breakthrough curves
The dispersivity is often estimated experimentally by assuming that the system is governed by Fickian Diffusion and fitting the exact solution ( 27) to a breakthrough curve, i.e. measurements of concentration over time at a fixed spatial location.In order to investigate the applicability and limitations of this approach when applied to a spatially variable dispersivity, here we compare our results for the 'true' solution to equations ( 18)-( 21) with the equivalent approximation assuming Fickian dispersion.
To do so, we can rewrite the Fickian solution (27) in terms of z and T .If the system has a dispersivity D = α 0 U x β but is instead assumed to have the Fickian dispersivity D F = α F U F , then the corresponding Fickian solution can be written as where is the ratio of dispersivities and U = U F /U is the ratio of average flow velocities.For a given value of β, we can fit (51) to the breakthrough curves from the numerical results described in Sect. 3 at a range of z locations.This fitting process can be done in two ways: either specifying the true average fluid velocity, U = 1, and determining σ alone, or fitting to the breakthrough curves with both U and σ as free, independent parameters.Once σ has been determined, fitting the scaling relation σ ∼ z γ (β) for each value of β then establishes a relationship between the inferred 'Fickian-fit' exponent γ discussed in Sect. 1 and the true exponent β.The two stages of this fitting process were carried out using MATLAB's nlinfit and fit routines, respectively.
Figure 5 shows breakthrough curves for values of β between 0 and 1.5, and the corresponding best-fit results from (51).The breakthrough curves are normalised by the late-time concentration value C ∞ (z) = C(z, T → ∞) at that z location -for β ≤ 1 this limiting value is just C ∞ = 1, but for β > 1 it is given by the steady solution (45) and so C ∞ < 1.The Fickian solution (51) is fit to these breakthrough curves in two cases: a '1-parameter' case with σ free and U = 1 fixed and a '2-parameter' case with both σ and U as free parameters.
Figure 5 shows that (51) accurately captures the β = 0 breakthrough curves, as expected.For β < 1 both the 1-parameter and 2-parameter cases produce good fits to the data, as the breakthrough curves are dominated by the late-time behaviour given by (38), which has a similar form to (51).For values of β close to 1 and above, however, the 1-parameter fit isn't capable of recreating the breakthrough curves we see herean apparent average flow velocity U F that varies with space is required.The best-fit velocity profiles U(z) from the 2-parameter approach are discussed in Appendix 1. Fig. 6 The results of the power-law fit for the dispersivity ratio σ .The coloured regions around each curve show their respective 95% confidence intervals.(a) The variation of the Fickian exponent γ with β.(b) The variation of the prefactor A, with the left y-axis corresponding to the 1-parameter result and the right y-axis to the 2-parameter result, for ease of comparison.(c) The normalised root-mean-square error (NRMSE) for the fitting to the normalised breakthrough curves, averaged over the number of time points and the number of breakthrough curves used.(d) The NRMSE for the power-law fitting to the best-fit σ values against z, averaged over the number of reference locations used and normalised by the mean of the σ values for each β We find that for small values of β the results are close to the true value, U = 1, but for larger values of β there is clear spatial dependence, with close-to-linear growth with z for β > 2.5.
For a given value of β, the diffusivity ratio σ (z) is well described by the power law σ ∼ Az γ , giving us a relation between the Fickian exponent γ and the true exponent β. Figure 6 shows the results of fitting this power law to σ for a range of β values, in both the 1-parameter and 2-parameter cases.The 95% confidence intervals are also plotted around each curve -these are larger for the 2-parameter case as there is now a trade-off between the values of σ and U when fitting, but the 2-parameter case still produces a better fit to the breakthrough curves as can be seen in Fig. 5 and Fig. 6(c).The prefactor A has differing behaviours in the two cases, but remains O(1).The Fickian exponent γ closely follows β for β < 1, due to the late-time behaviour (38) that mimics the β = 0 case.In the 1-parameter case, this continues for β > 1, while in the 2-parameter case after β ≈ 1.2, the corresponding value of γ begins to decrease but remains close to 1.As a result in the 2-parameter case, the relationship between β and γ is not one-to-one, and values of β greater than 1 all produce similar values of γ .For Fig. 6 each calculation for a particular value of β run until T = 5 , and breakthrough curves were then determined at spaced points between z = 10 0 and z = 10 3 .
To compare with results from experimental data, we express the Fickian dispersivity α F as α F ∼ α F,0 x γ .The curves in Fig. 6 provide a way to convert from the Table 1 A selection of reported values of α F,0 and γ , and the corresponding values of β and α 0 from the 1-parameter and 2-parameter curves in Fig. 6.In the 2-parameter case, there is the possibility of multiple(*) or zero( †) solutions, depending on the value of γ Source α F,0 γ β (1)  α (1) 0 β (2)   experimentally determined values of α F,0 and γ to A and β, and then to α 0 using the definitions of σ and z via A selection of reported α F,0 and γ values, and their corresponding α 0 and β values, is given in Table 1.
The results shown in Figs. 5 and 6 demonstrate that while this 'Fickian Fit' approach can be used as a suitable approximation in some cases, it fails to properly capture the full range of behaviour of a system with spatially dependent dispersivity such as this.The 1-parameter method estimates the value of β well, and by definition has the correct average fluid velocity of U = 1, but only produces the correct breakthrough curves small values of β.The method, conversely, reproduces all of the breakthrough curves well, at the expense of inaccurately estimating β and requiring unphysical velocity profiles for β > 1.To correctly determine the value of β from these breakthrough curves, the asymptotic solutions discussed in Sect. 4 should be used.

Discussion
The model discussed here presents a way to investigate the behaviour of a flow with power-law dispersivity, extending beyond the cases of β = 0 and β = 1 that have been considered by many before.The behaviours fall broadly into three classes: (i) the 'Fickian-like' case of 0 ≤ β < 1, where advection dominates at late time to give an error-function profile; (ii) the dividing case of β = 1, where advection and dispersion balance at all times; (iii) β > 1, where at late times dispersion is dominant and the concentration profile accumulates towards a steady solution that is significantly less than the inlet concentration away from the origin.
Wheatcraft and Tyler [10] demonstrated that a bundle of streamtubes in a onedimensional flow through a fractal medium of dimension d would have a dispersivity D ∼ x 2d−1 .Since for such a medium the fractal dimension will lie in the range 1 ≤ d ≤ 2, this corresponds to 1 ≤ β ≤ 3. Thus from a fractal point of view it is the β > 1 behaviour that is relevant for the one-dimensional flow, which has distinctly different behaviour to the Fickian case.Similarly, the case of radial flow can be transformed into one-dimensional flow, as shown in Sect.2.2, with the exponent δ from D ∼ r δ corresponding to β = 1 + δ/2.Thus Fickian radial flow corresponds to β = 1/2 one-dimensional flow, and β > 1 for δ > 0.
An often-used measure of the location of the advective front is where the concentration is half that at the inlet, i.e.C = 0.5.This is a reasonable assumption at late times for β = 0, as shown by (26) and Fig. 1(a), but is not necessarily the case for β > 0. Unlike in the Fickian case, dispersion is different on the two sides of the advective front for β > 0. This difference, characterised by the gradient of the dispersivity D (z), decreases as z → ∞ for β < 1 and so the concentration at the advective front C adv -never equal to 0.5 even for β = 0 -approaches 0.5 from below, as described by the late-time approximation (38).For β ≥ 1, however, D (z) increases as z → ∞ and C adv decreases from 0.5 at early times until it is controlled by the steady solution (45).These behaviours can be seen in the advective-frame plots in Fig. 1.For β > 1, the threshold of C = 0.5 would not be a good measure of the advective front, and for a location sufficiently far downstream would never in fact be reached.
The values of α 0 shown in Table 1 allow us to estimate the dimensional times that delineate the early-and late-time behaviours of the solutions shown in Sect. 3 by rewriting (17) as T .Based on the results shown in Fig. 1, we choose T = 10 −2 to mark the end of the early-time behaviour and T = 10 2 the start of the late-time behaviour, and take U = 1ms −1 as a representative value for the average flow velocity.We find that for β < 1 the values of α 1/(1−β) 0 are O (1) to O(10 1 ), or much smaller, so the early-time regime for a medium with β < 1 may not be of practical importance.In contrast for β > 1, the values of α 1/(1−β) 0 are O(10 1 ) to O(10 5 ), or much larger, so both the early-and late-time regimes may be physically relevant.
While this model extends a class of dispersivities that have previously been considered, it still has its limitations.It is unreasonable to assume that the dispersivity will continue to grow unbounded downstream of the injection point, as the properties of the porous medium will vary.While the exponent β can be related to the fractal dimension of the rock if this is known (which itself is only valid up to some cut-off length scale), we do not have a direct connection between it and the porosity and permeability of the porous medium.We can, however, extend our conclusions to a broader class of dispersivities D(z).If D (z) is a decreasing function then we would expect similar behaviour to the 'Fickian-like' case of 0 ≤ β < 1, with advection dominating at late times and a comparatively thin dispersive region around the advective front.On the other hand, if D (z) is an increasing function then we would expect behaviour like the β > 1 case, with the concentration profiles accumulating towards the steady solution with suitably chosen values of the constants A and B depending on the near-and far-field behaviour of D(z).
In this work, we have chosen to focus on a constant concentration inlet condition, C(0, T ) = 1, in order to describe the characteristics of a system such as this, with spatially dependent dispersivity, and to highlight the limitations of interpreting it as if it were Fickian.In practice, field and lab experiments typically inject tracers for a finite interval, introducing the timescale of tracer injection as a new parameter.The results presented here would then still be applicable for locations and times where the end of tracer injection has not yet been felt, but the inclusion of this more realistic inlet condition and this new parameter would be an important part of future work on this system.

Conclusions
In this paper, we have developed a mathematical model for anomalous solute transport in heterogeneous porous media with power-law dispersivity, D = α 0 U x β .This is motivated by experimental measurements of dispersivity at both the lab and field scale, along with the concept of dispersivity in fractal media as investigated by Wheatcraft and Tyler [10].This model includes the well-studied cases of Fickian dispersion (β = 0) and linear dispersivity (β = 1), extending them to the general case of β ≥ 0. We provide numerical results as well as exact analytical solutions or early-and late-time approximations.
These results, along with a simple scaling balance of the terms in the governing equation, demonstrate three classes of behaviour.For 0 ≤ β < 1 the system is 'Fickian-like', with dispersion dominating at early times and advection at late times.For β = 1 advection and dispersion remain in balance at all times, resulting in a self-similar solution.For β > 1 advection dominates at early times and dispersion at late times, and the concentration profiles tend to a non-trivial steady solution.This means that for β > 1 the concentration downstream will never reach that at the inlet.
In the Fickian case, at late times, the profile is close to symmetric about the advective front x = U t and the concentration there, C adv , is nearly 50% of that at the inlet.As a result of the spatial variability of the dispersivity for β > 0, increasing the degree of heterogeneity (i.e.β) causes the corresponding profile to become more asymmetric and the value of C adv to differ further from 50%; in fact for β > 1, C adv → 0 as time progresses.
Many experimental approaches have applied the Fickian solution to data and found that the inferred dispersivity scales with distance via an exponent γ .By following the same approach with our numerical results, we provide a direct connection between γ and the exponent in our model, β.
Here we have shown that, in general, solutions for the advection and dispersion of tracers in a porous medium with spatially correlated (or fractal) dispersivity lead to asymmetric concentration profiles.Such profiles, readily observable in both field and laboratory settings, yield important information on the nature and structure of the permeability field.It is our hope that the asymptotic solutions to the reduced model studied here may provide a computationally efficient way of classifying the degree of heterogeneity from measured concentration profiles, yielding new insights into longitudinal dispersion in aquifers.We recognise that this work is based on several simplifying assumptions and standardised conditions, for instance one-dimensional flow and a saturated medium, which may not necessarily be applicable to real scenarios.Therefore, in future work, a more concrete link to the underlying heterogeneities in rock properties needs to be included.Moreover, the current model is purely mathematical in nature and it would be helpful to see its experimental and field validation.7 Velocity profiles U (z) from the two-parameter Fickian Fitting process described in Sect. 5 for a range of values of β.Also shown are the best-fit power-law behaviours of U for small and large values of z, with their corresponding scaling exponents plotted in the lower right as a function of β velocity profiles are unphysical and a sign that while the 2-parameter Fickian Fit can accurately reproduce the breakthrough curves associated with this model, it is an incorrect description of the underlying physical system and its characteristic behaviour, as described in Sects.3 and 4.

Fig. 2 a
Fig. 2 a The dispersion-dominated approximate solution for a range of values of β < 1. b The advectiondominated approximate solution for all β ≥ 0

Fig. 4
Fig.4 The steady solution C s (z) for a range of values of β > 1

Fig. 5
Fig. 5 Breakthrough curves at a range of z locations for different values of β, normalised by the late-time concentration value C ∞ (z) at that z location.The breakthrough curves from this model are shown with solid curves.The best-fit Fickian curves (51) are also shown for σ free, U = 1 (dashed curve, •) and for both σ and U free (dotted curve, ×) Fig.7 Velocity profiles U (z) from the two-parameter Fickian Fitting process described in Sect. 5 for a range of values of β.Also shown are the best-fit power-law behaviours of U for small and large values of z, with their corresponding scaling exponents plotted in the lower right as a function of β