## Abstract

Solute transport in a fractured porous confined aquifer is modelled by using an equation with a fractional-in-time derivative of order *γ*, which may vary from 0 to 1. Accounting for non-Fickian diffusion into the surrounding rock mass, which is modelled by a fractional spatial derivative of order *α*, leads to the introduction of an additional fractional-in-time derivative of order *α*/(1+*α*) in the equation for solute transport. Closed-form solutions for solute concentrations in the aquifer and surrounding rocks are obtained for an arbitrary time-dependent source of contamination located at the inlet of the aquifer. Based on these solutions, different regimes of contaminant transport in aquifers with various physical properties are modelled and analysed.

## 1. Introduction

Aquifer contamination by contaminants of different nature (e.g. organic compounds, heavy metals, radioactive elements, etc.) is a widespread environmental problem. In many countries, it is common to use fractured bedrock aquifers as water supply and contamination of these aquifers is becoming a serious problem (Keller *et al*. 1995). The fractured porous aquifers are formed by porous rock matrixes of non-zero porosity dissected by a fractal-type network of fissures of high hydraulic conductivity (Fomin *et al*. 2003). For a fractured porous medium, fluid is stored in the porous elements and transported along the fissures. Water flow and solute transport by seeping groundwater are relatively slow, and it is not possible to carry out experiments over thousands of years and hundreds of meters of interest. Instead, one has to rely on models that describe the processes and mechanisms that would be dominant over long periods of time. It is therefore essential to understand the key processes well enough that credible predictions can be made using models based on well-established laws of nature. The difficulties in applying the basic laws of mass and momentum conservation arise mainly from the fact that the location, orientation and detailed hydraulic properties of the fractures cannot be measured in detail. The diffusion and absorption properties of the interior of the rock mass under natural stress also cannot be readily measured. The mixing processes of different water solutes in fractures and at intersections are not fully understood. This lack of measurement data makes it difficult to build models that account for all the processes. Simplifications are therefore made in an attempt to bring out the dominant processes.

Recently, Schumer *et al*. (2003) suggested a relatively simple mobile/immobile model with fractal retention times capable of simulating the anomalous character of solute concentration distributions for the flows in heterogeneous media of fractal geometry. In the field experiments carried out by Becker & Shapiro (2000), Haggerty *et al*. (2000) and Reimus *et al*. (2003) for the solute transport in highly heterogeneous media, the solute concentration profiles exhibited faster-than-Fickian growth rates, skewness and sharp leading edges. These effects cannot be predicted by the conventional mass transport equations. It was demonstrated in a number of publications (e.g. Meerschaert *et al*. 1999; Benson *et al*. 2000, 2001; Baeumer *et al*. 2001; Herrick *et al*. 2002; Schumer *et al*. 2003) that fractional differential equations can simulate the anomalous character of solute transport in highly heterogeneous media. Suggesting the inclusion of fractional-in-time derivative into the mass transport equation (in addition to the conventional derivative with respect to time), Schumer *et al*. (2003) referred to the conceptual model of multi-rate diffusion into immobile zones that had been described by Haggerty & Gorelick (1995), Cunningham *et al*. (1997), Carrera *et al*. (1998) and Haggerty *et al*. (2000).

In the present study, a mathematical model of contaminant transport in a confined, porous, fractured aquifer is derived and analysed. Through the use of fractional derivatives (in time and space), the model accounts for solute exchange between fissures and randomly distributed porous blocks in the aquifer as well as for the effects of contaminant transport retardation caused by non-Fickian diffusion into the porous blocks within the aquifer and into the confining rock.

For the case of an arbitrary time-dependent source of contamination located at the inlet of the aquifer, closed-form solutions for solute concentration in the aquifer and in the confining rock are obtained. Along with anomalous diffusion into the porous blocks within the aquifer, non-Fickian diffusion into the surrounding rock appears to be an important factor for retarding contamination. Interestingly, accounting for non-Fickian diffusion into the confining rock, which is modelled by a fractional-in-space derivative of order *α*, leads to emergence of an additional fractional-in-time derivative of order *α*/(1+*α*) in the solute transport equation. This implicitly justifies the applicability of the approach used by Schumer *et al*. (2003) for modelling diffusion in porous blocks and for describing the memory effects in solute transport by fractional-in-time derivatives. Hence, along with the non-fractured porous blocks within the aquifer, the rock that confines the aquifer from above and below can be treated as huge infinite porous blocks, which significantly contributes to the memory effects.

## 2. System model and analysis

A schematic of the fractured porous aquifer is presented in figure 1. Cartesian coordinates (*x*, *y*) are chosen in such a manner that fluid in the aquifer flows in the *x*-direction and that the coordinate *y* is directed upward. It can be readily shown that in the surrounding rocks confining the aquifer from below and above the gradient of solute concentration in the *x*-direction is much smaller than that in the direction orthogonal to the aquifer (Grisak & Pickens 1981; Tang *et al*. 1981; Tsang & Tsang 1987; Kennedy & Lennox 1995; Reimus *et al*. 2003; Rahman *et al*. 2004). Because the goal of the present study is to qualitatively estimate the effect of non-Fickian diffusion into the surrounding rocks on the solute transport in the aquifer (rather than to simulate the performance of concrete aquifer within the specific rock formation), it is assumed for simplicity that rocks below and above the aquifer have the same physical properties. The latter assumption makes the process symmetrical with regard to the median line of the aquifer *y*=−*h* (dashed line in figure 1). Therefore, the mass flux equals to 0 at *y*=−*h*, and the solution of the problem in the subdomain below this line is identical to the solution in the upper subdomain. Hence, only the upper half of the domain (*y*≥−*h*) can be considered. Thus, owing to the above assumptions, the differential equation describing the non-Fickian diffusive transport in the upper rocks can be presented in the following form:(2.1)where *D*_{1} is the effective diffusivity of porous media and the fractional derivative of order *α*, 0<*α*≤1, can be defined by means of Laplace transformation *L*, from the equation , which is equivalent to the Caputo definition (Samko *et al*. 1993)(2.2)The order of the fractional derivative in equation (2.1), *α*, is close to 0 for highly heterogeneous media and increases for more homogenous substances; *α*=1 for homogenous materials. In general, concentration *c*_{1} is a function of both spatial coordinates, *x* and *y*: *c*_{1}=*c*_{1}(*τ*, *x*, *y*). However, in equation (2.1), the derivative of *c*_{1} with respect to *x* is ignored. Dependence of *c*_{1} on *x* is a consequence of the boundary conditions on the rocks–aquifer interface (*y*=0), which couples *c*_{1} with the concentration in the aquifer, *c*_{2}=*c*_{2}(*τ*, *x*, *y*).

The advection–dispersion equation for the concentration of solute in a fractured porous aquifer, *c*_{2}, can be presented in the following form (Schumer *et al*. 2003):(2.3)where 0<*γ*, *λ*<1, *D*_{2} and are effective diffusivities, which account for dispersion and molecular diffusion in the aquifer in the *x*-direction and *y*-direction, respectively; *v* [*L*/*T*] is the fluid velocity; *τ* [*T*] is time; and *β* is the capacity coefficient that accounts for the ratio of immobile mass to mobile mass and for the physical properties of the blocks of porous rock within the aquifer (Haggerty & Gorelick 1995). The fractional-in-time, , and in-space, , , derivatives can be defined by means of Laplace transformation, *L*, as and , respectively, where *s* and *q* are the complex parameters of Laplace transformation or by the equivalent Caputo definition (see equation (2.2)). When 0<*γ*<1, the fractional-in-time derivative predicts the continuous transfer of contaminant from mobile to immobile phase (Schumer *et al*. 2003).

Integrating equation (2.3) over *y* and taking into account the mass flux continuity condition at *y*=0, as well as the condition of symmetry, at *y*=−*h*, and defining the mean value of concentration in the aquifer as , the advection–dispersion equation for mean concentration of solute in a fractured porous aquifer can be presented in the following form:(2.4)In order to convert equations (2.1) and (2.4) to non-dimensional form, proper characteristic scales must be defined. The scale for time, , represents the characteristic time for contaminant penetration in the rock matrix. The scale for the variable *x* along the aquifer, , is the characteristic distance of contaminant migration during the characteristic time-interval (0,*τ*_{m}). The scale for the *y*-coordinate is defined by the half-thickness of the aquifer, *h*. The initial concentration of solute at the inlet of the aquifer, *c*_{0}(0), where the source of contamination is located, can be used as the scale for solute concentration. Based on these scales, non-dimensional variables can be introduced as follows:(2.5)Substituting the non-dimensional variables from equation (2.5) into (2.4) yields the following:(2.6)It can be readily shown that, in many situations, contaminant transport within the aquifer in the direction of mainstream flow (*x*-direction) is characterized by relatively high values of the Peclet numbers *Pe* and, therefore, is predominantly determined by an advection mechanism (e.g. Neretnieks 1980; Grisak & Pickens 1981; Tang *et al*. 1981; Tsang & Tsang 1987; Haggerty & Gorelik 1995; Rahman *et al*. 2004). Hence, by imposing relevant boundary conditions, the non-dimensional mathematical model of the process can be presented in the following form:(2.7)(2.8)(2.9)(2.10)(2.11)(2.12)where *C*_{0}(*t*) is the non-dimensional concentration at the inlet of the aquifer. Note that assuming the concentration *C*_{1} on the rock–aquifer interface *Y*=0 to be equal to the mean concentration of solute in the aquifer, *C*, (boundary condition (2.11)) we slightly overestimate the value of concentration in the surrounding rocks. The boundary condition (2.11) is approximation of the more general mass transfer equation (Welty 2001): at *Y*=0, where the Sherwood number is greater than 10 for the regimes characterized by high Peclet numbers (owing to the high values of convective mass transfer coefficient *k*, which is much greater than the effective diffusivity of the rock matrix; Guo *et al*. 1999; McGuire *et al*. 2004).

Utilising Duhamel's theorem (Carslaw & Jaeger 1959), the solution for concentration, *C*_{1}, can be coupled with the concentration in the aquifer by the following equation:(2.13)where *u* is the solution of the auxiliary problem: it satisfies equation (2.7) with the initial condition (2.9), the boundary condition *u*=1 at *Y*=0 and vanishes at infinity (*Y*→∞, *u*→0).

Equation (2.13), in terms of mass fluxes, can be converted into the following form:(2.14)where and .

Applying the technique of Lie group analysis (Akcenov *et al*. 1994), it can be readily shown that the above formulated boundary-value problem, for function *u*, admits the following group of transformations:(2.15)where *ϵ* is the group parameter and *μ* is an arbitrary constant. The invariants of this group are(2.16)Therefore, the similarity solution in regard to the group of transformations (2.15) can be sought in an equation of the form of the equation *J*=*ω*(*η*), or, in other words, the solution *u* can be represented as(2.17)In order to satisfy the boundary condition *u*(0, *t*)=1, it is necessary to let *μ*=0 in equation (2.17). As a result, replacing *C*_{1} in equation (2.7) with *u*, (which is defined by formula (2.17), where *μ*=0), and accounting for the definition of fractional derivative, yields the following equation:(2.18)where *Γ*(*x*) is the Gamma function. The solution *ω* of equation (2.18) should satisfy the following boundary conditions:(2.19)(2.20)Applying Laplace transformation with respect to variable *η* to equation (2.18), and accounting for boundary condition (2.19), gives(2.21)where , and .

According to the definition of flux *Q*_{1}(*t*, *Y*) as the fractional derivative , and in view of formula (2.21), it can be readily seen that(2.22)or, in the particular case when *Y*=0,(2.23)From equation (2.21), it follows that(2.24)where *A* is a constant of integration.

Hence, accounting for the definition of the function *V*(2.25)Accounting for the well-documented properties of the Laplace transformation, the equations and , and boundary conditions (2.19) and (2.20), it can be readily shown that *A*=1 and that(2.26)Hence, as a result, formula (2.25) can be converted into the following form:(2.27)Accounting for the well-known relationship , where *Γ*(*x*, *y*) is the incomplete Gamma function, equations (2.27) and (2.26) can be rewritten in the following form:(2.28)(2.29)Equations (2.23) and (2.29) give(2.30)In order to find the solution for *ω*(*η*) of the boundary-value problem (2.18)–(2.20), it is convenient to first obtain the asymptotic forms of the function . Utilising the asymptotic formula for the incomplete Gamma function, when *s*→∞ (Abramowitz & Stegun 1972), yields(2.31)Calculating the inverse Laplace transformation, and accounting for formula (2.31), leads to the asymptotic representation of *ω*(*η*) for *η*<1,(2.32)The representation of *ω*(*η*) for *η*≫1 can be found analogously. Namely, asymptotically expanding the function and the incomplete Gamma function for *s*→0, in formula (2.28), yields(2.33)After the use of the inverse Laplace transformation, equation (2.33) leads to the following expression, which is valid for *η*≫1:(2.34)It is interesting to note that, although formula (2.32) is the asymptotic representation of the function *ω*(*η*) for small values of *η* (given that equation (2.31) was obtained for *s*→∞), the series that apparently follows from equation (2.32) for *M*→∞ will converge for all finite *η*. Therefore, formulae (2.32) and (2.34) define function *ω*(*η*) for all *η*∈[0,∞). Finally, accounting for formula (2.17) and using inverse substitution—that is, replacing *η* by in equations (2.32) and (2.34)—it can be concluded that the auxiliary solution *u* is obtained in closed form, and, consequently, the concentration of the solute in the confining rocks, *C*_{1}, can be evaluated from equation (2.13). The mass flux at *Y*=0 for this solution is defined by equations (2.14) and (2.30). In the particular case when diffusion is governed by Fick's law (*α*=1), the solution for *u* reduces to .

According to expression (2.30), equation (2.14), which couples the mass flux *Q*(*t*, 0) and *C*, can be presented in the following form:(2.35)Applying the Laplace transformation to formula (2.35) gives . Hence, in view of the definition of fractional derivative (Samko *et al*. 1993), equation (2.35) can be converted to the following form:(2.36)where, as can be readily seen, 0<*β*≤1/2 and *β*=1/2 for the particular case of Fickian diffusion (*α*=1).

Accounting for relationship (2.36), the boundary-value problem for solute transport in the aquifer, *C*, can be rewritten in the following form:(2.37)(2.38)(2.39)Note that in contrast to the model discussed in (Schumer *et al*. 2003), the problem formulated by equations (2.37)–(2.39) accounts for the variable-in-time source of contamination (boundary condition (2.39)) and contains an additional fractional-in-time derivative of order *β*=*α*/(*α*+1) (third term in the left-hand side of equation (2.37)), which models the non-Fickian diffusion into the surrounding porous rock.

As mentioned above, the mathematical model for the mean solute concentration (2.37)–(2.39) was derived assuming that the underlying rocks have the same physical properties as the upper rock strata. If this condition is not satisfied then the initial problem is not symmetrical relative to the line *y*=−*h*, then diffusion into the lower rocks should be also considered. In this case, equations (2.1) and (2.4) should be supplemented by the diffusion equation for the concentration *C*_{3} in the rocks below the aquifer. If the mass flux in the underlying rocks is determined by the fractional derivative of the same order as for the upper rocks, then it can be shown that accounting for diffusion in both upper and lower rocks with different effective diffusivities, *D*_{1} and *D*_{3}, leads to equation (2.37) with new constant factor (1+*D*_{3}/*D*_{1}) before the derivative .

The boundary-value problem (2.37)–(2.39) can be solved by using the Laplace transformation with respect to *t*. Using the Laplace transforms, equations (2.37)–(2.39) reduce to(2.40)(2.41)where .

The solution of equation (2.40), with the boundary condition (2.41), has the form(2.42)where(2.43)Calculating the inverse Laplace transform of formula (2.43) gives(2.44)From equations (2.42)–(2.44), calculating the inverse Laplace transform and accounting for Duhamel's theorem (Carslaw & Jaeger 1959), the solution of problem (2.37)–(2.39) can be presented as follows:(2.45)where *Χ*(*τ*) is a unit step function and function *φ* is defined by equation (2.44).

For *C*_{0}=1, which models a uniform source of contamination at *X*=0, solution (2.45) can be converted into the following form:(2.46)In the particular case where *γ*=1/2 and *β*=1/2, solution (2.46) obviously reduces to the well-documented formula (Neretnieks 1980; Tang *et al*. 1981)(2.47)The important practical case when contaminant leakage at *X*=0 occurs by accident and continues until a specific moment in time *t*^{*}, after which the source of contamination is eliminated and the flow of clean water resumes, can be modelled by assuming that *C*_{0}(*t*)=1 for *t* from the interval (0, *t*^{*}) and that *C*_{0}(*t*)=0 for *t*>*t*^{*}. Then, for this case, solution (2.45) can be converted to the form(2.48)If the solution for the concentration in the aquifer, *C*, is found, then the concentration in the surrounding rock, *C*_{1}, can be immediately obtained from equation (2.13), where *u* is defined by equations (2.17), (2.32) and (2.34). Solutions obtained for concentrations in the aquifer and surrounding rock matrix are rather simple and do not present any difficulties for numerical computations, especially, if such a universal and powerful tool as the computer algebra system Mathematica is applied. The accuracy of the numerical evaluations of the integrals in the solutions presented in general form can easily be verified by comparison to the simple solution (2.47) obtained for the particular case when *C*_{0}=1 and *β*=*γ*=1/2.

## 3. Results and discussion

Prior to numerical computations, some conclusions can be made just on the basis of the governing equations derived above. For instance, the additional fractional-in-time derivative of order *β*=*α*/(*α*+1) in equation (2.37), which appeared as a result of accounting for non-Fickian diffusion in the surrounding rock (see equation (2.7)), also substantiates the approach of Schumer *et al*. (2003). According to this approach, the interaction of mobile and immobile phases in the aquifer can be modelled by the fractional-in-time derivative (the second term in the left-hand side of equations (2.8) and (2.37)), in which the capacity parameter, *b*, is determined by the ratio of the masses of the immobile and mobile phases at equilibrium and by the diffusivity, porosity, geometry and size of the blocks, and in which the order of derivative *γ* is determined by the peculiarities of the distribution of pores and micro-cracks in the blocks, by the adopted model of mass transfer through the blocks and by the distribution of these blocks in the aquifer (Haggerty & Gorelick 1995).

It is also worth noting that the order *β* of the third derivative of time in the left-hand side of equation (2.37), which accounts for diffusive transport into the confining rock, reaches its maximum value, 1/2, in case of Fickian diffusion (*α*=1), and for non-Fickian diffusion *β*<1/2. Extending the latter observation to the fractional-in-time derivative of order *γ*, which accounts for mass transfer into the immobile phase (porous blocks within the aquifer), leads to the conclusion that if mass transport into the porous blocks within the aquifer is determined only by diffusion, then the order of the fractional derivative *γ* should not exceed 1/2. Of course, the latter conclusion is valid only for relatively short periods of time of contamination and only if the sizes of the porous blocks within the fractured aquifer are relatively large, so that the contamination front from one side of the porous block in the aquifer does not meet the contamination front moving in from the other side of the block (e.g. if the concentration fields in the blocks are disturbed in the narrow region in the vicinity of their boundaries).

From the solutions obtained for *C* and *C*_{1} (equations (2.44), (2.45), (2.13) and (2.17)), it can be readily seen that for a uniform source of contamination at the inlet of the aquifer, *C*_{0}(*t*)=1, the long-time concentrations in the aquifer and surrounding rock (for finite *X* and *Y*) behave like the following:(3.1)(3.2)Hence, if contamination continues over a long time and the source of contamination is not eventually eliminated, then concentrations in the aquifer and in the neighbouring rock masses tend to match the concentration at the inlet. For *γ*<*β*, the distribution of concentration is determined by diffusion in the porous blocks in the aquifer, otherwise (for *γ*>*β*) the effects of diffusion in the surrounding rock will dominate. In another particular case of practical importance, contamination may occur due to the instantaneous injection of fluid at *X*=0. This can be modelled by a Dirac *δ*-function assuming that *C*_{0}=*δ*(*t*). In this case, the long-time asymptotic solutions (for finite *X* and *Y*) decay like(3.3)(3.4)Again, in this case, for *γ*<*β*, the decay of concentration takes place as a result of diffusion through the porous blocks within the aquifer, whereas for *γ*>*β*, it happens mostly due to diffusion into the surrounding rock. Note, however, that in all situations, equations (3.3) and (3.4) predict a more rapid asymptotic decline of concentration than the long-time decline of concentration in the case of constant inlet contamination described by equations (3.1) and (3.2). It is also interesting to note that if the effect of diffusion in the surrounding rock is ignored (i.e. if the first term in the right-hand-side of equation (3.3) is omitted), then equation (3.3) coincides with the asymptotic representation for *C* obtained by Schumer *et al*. (2003) for non-confined flow.

Figure 2 illustrates the distribution of concentration *C*_{1} in the confining rock matrix with respect to coordinate *Y* (perpendicular to the aquifer) from a constant source of contamination, *C*_{0}=1, located at the inlet of the aquifer, *X*=0, for various durations of contamination period. As it can be seen, near the aquifer, in the case of non-Fickian diffusive transport (*α*=0.5), for all periods of time (*t*=10; 50; and 100) the corresponding concentrations decline faster than in the case of Fick's model. For longer distances, the contaminant concentration for non-Fickian transport is significantly higher than can be predicted by Fick's law. In other words, for long periods of time, smaller values of *α* simulate solute diffusion into the surrounding rock in greater volumes. This leads to greater retardation of contaminant transport in the aquifer in the streamwise direction. On the other hand, deeper penetration of the contaminant by diffusion into the surrounding rock matrix may, over time, lead to contamination of the upper subterranean strata and may eventually affect areas in the vicinity of the Earth's surface.

Distributions of concentrations in the aquifer, with respect to the stream-wise distance *X* from a constant source of contamination, *C*_{0}=1 at *X*=0, for various orders of derivatives *γ* and *α* in equation (2.33), and for a contaminating period duration of *t*=15, were computed by equation (2.46) and are shown in figure 3*a*–*d*. For relatively long periods of time (e.g. for *t*=15), the corresponding concentration varies according the asymptotic formula (3.1). Namely, for example, at *α*=1 (i.e. *β*=0.5), the last term in the right-hand side of expression (3.1) for *γ*=0.2 is greater than for *γ*=0.5 and therefore the concentrations decreased, whereas for *γ*=0.8, this term becomes smaller, which leads to a higher concentration (see also figure 3*a*). In other words, smaller values of *γ* simulate solute transition into the immobile phase in greater volumes; that is penetration of greater volumes of solute into the porous blocks within the aquifer. A similar effect can be observed by varying *α* when *γ* is fixed (see figure 3*b*,*c*, where *γ*=0.5). In this case, smaller values of *α* model a more intensive diffusion into the confining rock, and, therefore, the concentration of solute in the mobile phase decreases. The results presented in figure 3*d* indicate that reduction of parameter *b* leads to less transition of solute into the immobile phase in the aquifer, which results in a higher concentration, *C*, of solute in the mobile phase. In all cases the highest concentration of contaminant is observed when diffusion in the confining rock is modelled by Fick's law (*α*=1) and transition into the immobile phase is ignored (*b*=0), which is illustrated by the dashed curves in figure 3*a–d*.

The break-through curves presented in figures 4*a*,*b* and 5*a*,*b* illustrate the variation of concentration, *C*, in the aquifer over time at point *X*=3. These results are computed by using the close-form solution (2.48) and (2.44), which models two regimes: (i) the regime of contamination from a source, located at *X*=0, over a finite period of time (0, *t*^{*}) and (ii) the consequent regime of clean water flow that begins at *t*=*t*^{*}. Comparison of the corresponding plots in figure 4*a* shows that smaller values of parameter *γ* lead to lower concentrations in the aquifer. These differences in concentration initially become more pronounced over time, reach their maxima at certain times and become less pronounced over longer times. As mentioned above, smaller concentrations are typical for smaller values of *γ*, and this can be attributed, for example, to more efficient diffusion of the contaminant into the porous blocks. A similar effect can be observed in figure 4*b*, which illustrate the case where diffusion in the blocks is ignored (*b*=0), and the retardation of contaminant transport is determined by the diffusion into the surrounding rock. Here also, anomalous diffusion (i.e. small values of *α*) leads to stronger penetration of the contaminant into the surrounding rock and, as a consequence, to a lower concentration in the aquifer. A similar behaviour can be seen in figure 5*a*, in which the solute transition into the immobile phase (diffusion into the porous blocks within the aquifer) is not ignored (*b*=0.5). Results presented in figure 5*b* for different values of parameter *b* indicate that reducing this parameter leads to a higher concentration of solute in the aquifer. A smaller *b* corresponds to a smaller fraction of immobile phase; that is, a smaller number of porous blocks or a lower porosity for these blocks, which reduces the effective diffusivity of the blocks within the aquifer.

The results presented in figure 6*a*,*b* illustrate lateral distributions of average aquifer solute concentration *C* for different times. Solid lines correspond to the stage of continuing contamination and dashed lines to the regime of clean water flow after elimination of the source of contamination. In the aquifer, similarly to what has been observed for the concentration in the confining rocks (see figure 2), the higher concentrations *C* are obtained for the higher orders *γ* of time derivative (therefore, *γ*=0.8 corresponds to the less intensive diffusion into the porous blocks within the aquifer). Smaller values of *γ* simulate diffusion into the blocks, which form the aquifer, in greater volumes (greater transition of the contaminant into an immobile phase), leading to the effect of greater retardation of contaminant transport in the aquifer.

## 4. Conclusions

Using the approach suggested by Schumer *et al*. (2003), the anomalous mass transport in a fractured porous aquifer and surrounding rock is modelled, and closed-form solutions for the corresponding governing equations with fractional derivatives (in time and in space) are obtained.

The effect of non-Fickian diffusion in the surrounding rock mass can be modelled by including a fractional-in-time derivative of order *β*=*α*/(*α*+1) (where 0≤*α*≤1 is the order of a derivative that models non-Fickian anomalous diffusion within the rock matrix) in the equation for solute transport in the aquifer. This result once again confirms the applicability of the approach used by Schumer *et al*. (2003) for modelling solute transport in a fractured porous aquifer by equations with fractional-in-time derivatives.

Introducing a fractional-in-time derivative of arbitrary order into the equation of mass transport for a fractured porous medium leads to an efficient model for describing the mechanisms of fluid–rock interaction within the fractured aquifer of a complex geological structure.

The parameters that simulate different regimes of fluid–rock interaction and different properties of porous blocks within the aquifer are indicated, and their effects are analysed. For example, the long-time behaviour of concentration distribution in the aquifer and in the rock matrix suggests that low values of *γ* (with the same capacity coefficient *b*) lead to decreased concentrations in both the aquifer and in the surrounding rock, whereas lower values of coefficient *b* result in higher solute concentrations.

## Footnotes

- Received July 13, 2004.
- Accepted April 1, 2005.

- © 2005 The Royal Society