Fractal advection-dispersion equation for groundwater transport in fractured aquifers with self-similarities

Groundwater transport within a fractured aquifer with a fractal nature exhibiting self-similarity cannot be accurately simulated by the classical Fickian advection-dispersion transport equation without a detailed characterisation of the fracture network and heterogeneity of the system. Because the information to characterise such a system to an appropriate level of detail is often not available, most applications fail to accurately simulate the observed contaminant transport. In response to the current limitations of transport modelling using the advection-dispersion equation, especially in fractured media, a fractal advection-dispersion groundwater transport equation is developed. Fractal differentiation is discussed in terms of the fractal derivative and the fractal integral. The fractal derivative is commonly applied and known, yet the fractal integral is developed in this paper along with the appropriate theorem and proof. The numerical approximation of the fractal derivative and integral are given, where Simpson’s 3/8Rule and Boole’s Rule for numerical integration are applied for the fractal integral, and the forward and Crank-Nicolson finite difference schemes are applied for the fractal derivative. Upon the given foundation, the fractal advection-dispersion transport equation is formulated to develop a new groundwater transport model. The qualitative properties of the fractal advection-dispersion equation are investigated to determine boundedness, existence and uniqueness of the solution. To validate the developed fractal advection-dispersion equation, a numerical simulation is performed with different fractal dimensions to demonstrate the applicability of the model to fractal groundwater systems. The numerical simulations found that anomalous diffusion could be better modelled by incorporating a fractal dimension, especially where limited information is available on preferential pathways.


Introduction
The classical Fickian advection-dispersion transport equation application to groundwater systems to model transport relies on the availability of information to characterise the physical system in its entirety. However, where preferential flow pathways, such as fractures, faults or dykes, are present there is often insufficient information to accurately characterise the heterogeneity of these systems. In these environments with an incomplete characterisation of heterogeneity, the classical Fickian advection-dispersion transport equation is unable to accurately simulate the observed plume [1,2]. Characterising the heterogeneity of a system is further complicated by the fact that dispersivity increases with the scale of measurement [1][2][3][4][5]. Traditional solutions to the problem of dispersivity scale-dependency consist in establishing scaling relationships [2], or alternatively models are used to incorporate fractures (or preferential flow pathways), namely the dual porosity model; random array of fractures; multiple continuum; and others [6,7]. However, these models cannot account for the fractal nature of fractured systems.
Fractal geometry and fractured systems were originally related in 1985 for a nuclear waste disposal site investigation [8], and a number of authors have corroborated this relation [7,[9][10][11]. Fractal geometry is based on the principle 2 Fractal differentiation

Fractal derivative
The fractal derivative differs from the conventional integer-order derivative, where the integer-order derivative represents the change of a function (dependent variable) with the change of another quantity (independent variable) in ordinary space, while the fractal derivative represents the ratio of change of two quantities in fractal space. The fractal derivative can be defined by transforming the standard integer dimensional space-time (u, t) into fractal time (u, t α ) [14,16,17]: where α denotes fractal dimension of time.
Alternatively, the fractal derivative can be defined by transforming standard integer dimensional space-time (u, t) into fractal space-time (u β , t α ) [14]: where β denotes the fractal dimension of space.

Fractal integral
The fractal integral or antiderivative can be determined by considering the fractal time derivative of a function, and assuming that the derivative is known and denoted by a function u(t): where α denotes fractal dimension of time.
To further solve eq. (5) the Laplace transform (L) is applied Remembering the Laplace transform of an integer-order derivative: Rearranging, simplifying and dividing by the Laplace space (s): The inverse Laplace transform is applied to convert back from the Laplace domain: The convolution theorem becomes useful for calculations, where the convolution theorem gives the inverse transform of the product of two transforms. Considering two functions f and g, which are piecewise continuous on t ≥ 0, the Laplace transform is [18]: The convolution of f and g is (11) and, the inverse Laplace transform is thus: The Laplace transform is additive, but not multiplicative. This means that the Laplace transform of a product is not the product of the Laplace transforms. The convolution theorem gives the transform required to get a product of Laplace transforms, i.e. the convolution [18]. Applying the convolution theorem to eq. (9): Thus, the fractal integral can be defined as Theorem 1 [19]: Let the function (f ) be differentiable in an open interval I. Then the fractal integral of the function (f ) is given as . Consider the definition of the fractal integral (eq. (15)): Remembering that the fractal derivative can be expressed as d (3) and eq. (4)): Simplifying, Applying the integral from t to 0: Secondly, the following is considered for the proof: Let the new function F be defined as Applying the definition of a fractal derivative (eq. (1)): Remembering that the fractal derivative can be expressed as d (3) and eq. (4)): Consider the definition of the fractal integral (eq. (15)): Simplifying using the fact that d dt ( This completes the proof.

New groundwater transport model: fractal advection-dispersion equation
The traditional advection-dispersion equation is used to describe the transport of non-reactive contaminants or tracers, which is founded on an analogy to Fick's law of diffusion. Non-Fickian transport or diffusion, also termed anomalous diffusion, is any transport which is not adequately described by the traditional advection-dispersion equation. Anomalous transport can be defined by non-Gaussian leading or trailing tails of a plume or break-through curve from a point source, or nonlinear growth of the centred second moment. Subdivisions of anomalous diffusion include superdiffusion and subdiffusion, where superdiffusion is characterised by a growth rate faster than linear growth (i.e. faster movement found in preferential flow pathways), and subdiffusion is characterised by a growth rate slower than linear growth (i.e. slower movement in low-permeability zones) [20][21][22][23]. These shortcomings of traditional advection-dispersion equation have led to alternative non-local conceptualisations of flow and transport, and numerous methods for addressing scale and space-time dependencies. These alternative methods include stochastic averaging of the classical advection-dispersion equation, the multiple-rate mass transfer method, the continuous time random walk method, the time fractional advection-dispersion equation method, the space fractional advection-dispersion equation method, and others [20]. In these alternative methods, the dispersive state is allowed to vary between superdiffusion, subdiffusion and normal diffusion, referred to as transient dispersion [22]. However, each method might be formulated for a specific transition and thus not appropriate for all types of transient dispersion. The fractional advection-dispersion equation formulations have proven successful in describing non-Fickian transport, yet for three-dimensional applications a scale-dependent dispersivity problem has been found [15,24], similar to the traditional advection-dispersion equation.

Motivation
Groundwater transport within a fractured aquifer with a fractal nature exhibiting self-similarity cannot be accurately simulated by the classical Fickian advection-dispersion transport equation without a detailed characterisation of the fracture network and heterogeneity of the system. Because the information to characterise such a system to an appropriate level of detail is often not available, most applications fail to accurately simulate the observed contaminant transport. Two approaches have been used to overcome this limitation of the classical Fickian advection-dispersion transport equation, one being methods to improve characterisation of the physical system, and secondly alternative mathematical formulations such as the fractional derivative. Yet, neither of these approaches can account for the fractal nature of fractured systems.
Developments in the theory of fractal geometry have resulted in numerous applications, to many branches of science, especially to problems stemming from unstable phenomena [6,[25][26][27][28][29]. Fractal geometry incorporates the complexity of the natural system by describing an object by the repetition of a given algorithmic rule or pattern over a multitude of separate length scales. Fractured-rock aquifers have proven to be an appropriate candidate for a fractal description [6], and the ability of a fractal geometry to describe complexity on different scales potentially provides the ability to describe scale-dependent dispersivity.
In response to the current limitations of transport modelling using the advection-dispersion equation, especially in fractured media, a fractal advection-dispersion groundwater transport equation is developed.

Fractal formulation of the advection-dispersion transport equation
To include the effects of self-similarity inherent in groundwater transport with preferential pathways formed by fractures, the mathematical formulation of the one-dimensional advection-dispersion equation is altered to include a fractal space component: where α denotes the fractal dimension To express the fractal ADE in terms of integer-order dimensions, the defined property of the fractal derivative Without the loss of generality, the advection transport term is considered Then, the dispersion transport term is considered Substituting the advection and dispersion terms back into eq. (26): Thus, eq. (31) simplifies to where V α F (x) is the velocity with fractal dimension with respect to x, and D α F (x) is the dispersion coefficient with fractal dimension with respect to x.
Remark: The fractal advection-dispersion equation should revert back to the classical formulation when α = 1.
Considering the fractal dimension velocity term (α = 1): Considering the fractal dimension dispersion term (α = 1): When α = 1, the fractal advection-dispersion equation reverts back to the classical formulation of the advectiondispersion equation.

Qualitative properties
The boundedness, existence and uniqueness of the fractal transport model (eq. (32)) is first determined before moving to the numerical approximation in the following section. The Picard-Lindelöf theorem will be used to achieve this.

Lipschitz condition boundedness for partial differential
The Lipschitz condition is used for a bound on the modulus of continuity for a function. A function f satisfies a Lipschitz condition on a set B, if there is a constant M ≥ 0, such that A function which has a bounded first derivative is considered Lipschitz, and this has been assumed to be true for partial differential operators. A theorem and proof are proposed here for the partial differential equations.
Since f, g ∈ B, and similarly f − g ∈ B, then by assumption a positive constant M can be found, such that And, multiplying by f −g This completes the proof.

Fixed-point theorem for existence and uniqueness
An integral is applied to both sides of the fractal transport equation (eq. (32)), to obtain Let a new function F C(x, t) be expressed as and, state the Lipschitz condition eq. (41): Without the loss of generality, consider the term left of the inequality sign, and substitute the function F C(x, t): Simplifying and factorising Applying the triangular inequality and remembering Applying the proven theorem in sect. 3.3.1, Thus, the fractal transport equation does uphold the Lipschitz condition.
To facilitate the evaluation of the fractal transport equation, a new function is introduced Integrating on both sides of eq. (51), A set is created in which to evaluate the fractal transport equation where The Banach fixed-point theorem is applied by introducing the norm of the supremum (statistical limit of a set) for Remembering the physical meaning of the fractal advection dispersion equation and the nature of groundwater transport, the spatial component (b) of an aquifer domain is such that the simulated plume when the initial concentration source is removed cannot be greater than the aquifer total extent: Remembering eq. (50), Taking the integral out of the normalised term Considering the individual components Considering the physical nature of the components in eq. (59), the maximum value of the fractal velocity and dispersivity terms can be used to define boundedness, and the nature of a derivative implies that it is bounded and can be represented by Φ: Now, the terms can be moved out from the integral Taking the integral of the change in continuous time (dτ ), and incorporating the maximum time for which the contaminant could be found, Remembering eq. (56) as the condition for a well-defined problem, The operator is thus well defined, where Furthermore, following the equality in eq. (49), the contraction can be defined as (K < 1): Thus, the defined operator is well defined and is a contraction for the condition: Thus, according to the Banach fixed-point theorem, the developed fractal transport advection dispersion equation has a unique solution.

Numerical approximation of the new fractal transport model
Numerical approximation methods for the newly developed fractal advection-dispersion equation are presented, where Simpson's 3/8 Rule and Boole's Rule for numerical integration are applied for the fractal integral, and the forward and Crank-Nicolson finite difference schemes are applied for the fractal derivative.

Fractal derivative: finite difference numerical approximation
The fractal derivative is a local operator which allows for the application of traditional finite difference approximation methods for time and space. The fractal advection-dispersion equation is numerically approximated using the traditional forward finite difference and the Crank-Nicolson finite difference methods.

Forward finite difference scheme
Explicit forward differences in time and the second-order derivative is applied to eq. (32): where i denotes a one-dimensional grid-centered framework with conventional unit vector (i), V α F is the fractal dimension velocity term, and D α F is the fractal dimension dispersion term. Rearranging c k+1 Simplifying by using constants a 1 , b 1 , c 1 , and d 1 where A finite difference scheme is considered stable if the errors incurred at a discrete time step are not propagated throughout the simulation. A Von Neumann stability analysis is conducted to check the stability of finite difference schemes applied to the fractal advection-dispersion equation. The error incurred in the numerical approximation can be defined as where λ k i is the approximation or round-off error, N k i is the numerical solution, and c k i is the exact solution. Considering eq. (69) in terms of the recurrence relation for the error: The error for each discrete point in time and space has been determined as Simplifying and rearranging The amplification factor (G) can be defined as where, for the solution to remain bounded the condition |G| ≤ 1. Thus eq. (78) expresses the stability criteria for the fractal advection-dispersion equation, where |e aΔt | ≤ 1, and remembering that a 1 = 1 Δt : Applying known identities and the absolute value of a complex number: Equation (81) is an expression of the stability criteria for eq. (67).

Crank-Nicolson finite difference scheme
When the Crank-Nicolson scheme is applied Rearranging, Simplifying by using constants a 1 , b 1 , c 1 , and d 1 where A Von Neumann stability analysis is also conducted to check the stability of the finite difference approximation applied to eq. (84), remembering eqs. (70) to (75), and additionally: To obtain Simplifying and rearranging Applying known identities, simplifying and applying the absolute value of a complex number, where z = a + ib and |z| = √ a 2 + b 2 : where Simplifying, .
Equation (92) is an expression of the stability criteria for eq. (82).

Fractal integral: numerical integration approximation
The fractal derivative is a local operator, and thus the numerical solution of fractal derivative equations can be approximated using standard numerical techniques for the integer-order derivative equations [14]. For the numerical integration of the fractal integral, standard integer-order techniques are also utilized, yet instead of the traditionally used Trapezoid Rule for numerical integration the more rigorous Simpson's 3/8 Rule and Boole's Rule are applied.
Equation (100) is the numerical approximation of the fractal integral using Simpson's 3/8 Rule for numerical integration.

Numerical simulation for different fractal dimensions
The developed fractal advection-dispersion equation is numerically simulated for a generic groundwater transport problem to investigate the effects of varying the fractal dimension. The one-dimensional first-type boundary condition transport problem aims to determine the concentration of a contaminant at a specific distance from the source. The concentration of the contaminant (C 0 ) is set at 1 000 mg/l, and defined as a continuous source at x = 0: A uniform groundwater flow velocity and longitudinal dispersivity are assigned to the receiving aquifer. The groundwater velocity was set to 0.5 m/d and the longitudinal dispersivity was set to 0.3 m. The initial concentration in the aquifer is set to 0 mg/l and thus the initial conditions are defined as The contaminant movement is set to have a finite distance and the boundary condition is thus defined as: The solution of this generic transport problem is simulated using one of the four numerical solution methods described in this paper for the developed fractal advection-dispersion equation, namely the Crank-Nicolson finite difference approximation method (eq. (82)-eq. (84)). The fractal advection-dispersion equation simplifies to the classical advectiondispersion equation when α = 1. The finite difference numerical solution of the classical advection-dispersion equation (α = 1) is displayed in fig. 1, which will serve as a base for comparison to the fractal solutions (0 < α < 1). The solution of the classical advection-dispersion equation for the defined transport problem illustrates how the continuous contaminant source is periodically transported under the constant velocity flowing in the positive x-direction. At a time of 100 days, the simulated concentration at a distance of 400 m from the source is 0 mg/l, and after 1 000 days the concentration at the same distance is approximately 600 mg/l. The effect of varying the fractal dimension (α) on the solution to the fractal advection-dispersion equation is explored using the Crank-Nicolson finite difference approximation method. The fractal dimension (α) is varied between the limits (0 < α < 1) with an increment of 0.1. The one-dimensional solution for each of these realisations of the fractal dimension is illustrated in figs. 2 and 3.  Adding the fractal dimension to the simulation of the groundwater transport in the defined one-dimensional transport problem, has the effect of transporting the contaminant preferentially along the direction of groundwater flow. For example, the simulated concentration at a distance of 400 m from the sources is approximately 600 mg/l at the initial time already, compared to the classical simulation where this concentration was only reached at the end of the simulation period of 1 000 days. This transport could be characterised as superdiffusion, which is characterised by a growth rate faster than linear growth (i.e. faster movement found in preferential flow pathways). Interestingly, the fractal dimension additionally simulates transport in a direction parallel to the defined groundwater flow direction, with clear self-similarity. This additional preferential pathway can be seen from fractal dimensions of 0.2 to 0.9, with the definition of these pathways becoming clearer with an increasing fractal dimension number. Considering the solution with a fractal dimension of α = 0.5, there are clear preferential zones of groundwater transport.
It is important to note that a constant unidirectional groundwater velocity was assigned to the model, yet preferential pathways parallel to the assigned groundwater velocity direction are realised with the addition of the fractal dimension. From this result, it can be hypothesised that the fractal formulation of the advection-dispersion equation could simulate preferential pathways implicitly. Thus, simulate the effect of fractures or faults without explicitly including the location and specific details of the fracture or fault system in the model. In this way, contaminant transport showing signs of anomalous diffusion, could be better modelled by incorporating a fractal dimension. Especially, where limited information is available on the preferential pathway causing the discrepancy.

Conclusions
It has been discussed that groundwater transport within a fractured aquifer with a fractal nature cannot be accurately simulated by the classical Fickian advection-dispersion transport equation without a detailed characterisation of the fracture network and heterogeneity of the system. The use of the classical advection-dispersion equation tends to inaccurately simulate the observed contaminant transport because the information to characterise a fractured system to an appropriate level of detail is often not available. This is due to the fact that heterogeneity is not explicitly incorporated into the classical advection-dispersion equation, but rather a constant velocity and dispersivity is considered at each point and heterogeneity is integrated externally. In response to the current limitations of transport modelling using the advection-dispersion equation, especially in fractured media, a fractal advection-dispersion groundwater transport equation was developed and four methods to solve the fractal advection-dispersion equation are described, namely forward finite differences and Crank-Nicolson finite differences for the fractal derivative equation, and Simpson's 3/8 and Boole's numerical integration for the fractal integral equation.
The developed fractal advection-dispersion equation was numerically simulated for a generic groundwater transport model to investigate the effects of varying the fractal dimension in a one-dimensional groundwater flow system. The simulations provided evidence that the fractal formulation of the advection-dispersion equation could possibly model superdiffusion, which is characterised by a growth rate faster than linear growth (i.e. faster movement found in preferential flow pathways). Additionally, the fractal dimension simulates transport in a direction parallel to the defined groundwater flow direction, with self-similarity. It was thus theorised that the fractal formulation of the advection-dispersion equation could simulate preferential pathways implicitly. Where the effect of fractures or faults could be modelled without explicitly including the location and specific details of the fracture or fault system in the model. In this way, anomalous transport could be better modelled by incorporating a fractal dimension, especially, where limited information is available on the preferential pathway causing the discrepancy.
Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.