# Impact of Hydraulic Conductivity on Solute Transport in Highly Heterogeneous Aquifer

Kaili Wang, Guanhua Huang

► To cite this version:

Kaili Wang, Guanhua Huang. Impact of Hydraulic Conductivity on Solute Transport in Highly Heterogeneous Aquifer. 4th Conference on Computer and Computing Technologies in Agriculture (CCTA), Oct 2010, Nanchang, China. pp.643-655, 10.1007/978-3-642-18333-1_78. hal-01559611

HAL Id: hal-01559611

https://inria.hal.science/hal-01559611

Submitted on 10 Jul 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Impact of Permeability on Flow and Solute Transport in Highly Heterogeneous Aquifer

Kaili Wang¹,² Guanhua Huang¹,²*

1. College of Water Conservancy and Civil Engineering, China Agricultural University, Beijing, 100083, P.R. China

2. Chinese-Israeli International Center for Research and Training in Agriculture, China Agricultural University, Beijing, 100083, P. R. China

**Abstract.** The impact of permeability on solute transport process in a highly heterogeneous aquifer was analyzed using the Monte-Carlo method. The logarithm of the permeability (lnK) of aquifer was considered as a non-stationary field with increments being a truncated fractional Lévy motion (fLm) generated using the SRA3DC code. MODFLOW and MT3DMS code were used to solve the flow and solute equations, respectively. Results indicate that larger C lead to a more heterogeneous permeability field. As C increases, solute plumes show more significant anomaly with sharper leading edge and wider tailing edge. The solute plume extends and its second spatial moments increase as C increases, while the first spatial moments of the solute plume are independent of C values. The longitudinal macrodispersivity is scale-dependent and increases as a power law function of time. Increasing C results in an increase in longitudinal macrodispersivity. It was concluded that the nature of transport process is highly dependent on the heterogeneity of the permeability field.

**Keywords:** permeability, transport, Monte-Carlo, heterogeneous, moments

## 1. Introduction

Groundwater contamination has become one of the most important environmental issues all over the world. It is necessary to predict flow and contaminant spreading in

* Corresponding author: Guanhua Huang, E-mail: 53422446@163.com
the subsurface for the control of groundwater quality. However, the heterogeneity of porous media and incomplete knowledge of data information lead to difficulty in the estimation of hydraulic properties and geophysical variables, and thus, bring about difficulty in estimating or predicting subsurface flow and transport. The stochastic methods have been developed and used to deal with these difficulties [1]. There are mainly two categories of stochastic methods in widespread use including the moment equation method (MEM) and the Monte-Carlo method (MCM). Compared with the MEM, the MCM is conceptually and computationally straightforward, and most importantly, once the properties' distributions are established in details, the MCM can be applied to highly heterogeneous media. Up to present, the MCM has been widely used in researches on flow and transport process with the more accurate subsurface simulation models and the availability of high-speed supercomputers [2-6].

In the MCM simulation, the first step is to determine the distributions of the properties of geological formations. Researches have shown that the permeability (K) has a relative large spatial variability with its maximum value several orders being larger than its minimum value, while the spatial variability of porosity and other geophysical parameters are relatively small [7]. In classical stochastic theories, K is treated as a lognormal, stationary and spatially correlated random function, which can be characterized by its mean, variance and single finite-range integral scale [1]. It has been proved that such a treatment can be only used for those media with relatively small heterogeneities [8]. In recent years, several researches have shown that the geological media can be highly heterogeneous and exhibit very large and abrupt changes in K field [9, 10]. And such a hydraulic field with continuously evolving heterogeneities can be described using the fractal scaling model [4, 11, and 12]. One of the commonly used fractal model to characterize the non-stationary K field is the fractional Brownian motion (fBm), with which the probability density function (PDF) of lnK increments for any given lag are assumed to be Gaussian. However, researches show that the lnK increments of heterogeneous formations are often non-Gaussian with an increased peak value around its mean and heavy (slow-decaying) tails [10]. These non-Gaussian behaviors can be characterized using the probability distributions of the Lévy-stable family, which is referred to as the fractional Lévy motion (fLm). Compared with the Gaussian distribution, the Lévy distribution can better mimic the sharp property contrasts associated to the geological formations. However, the heavy tails also give the Lévy-stable distribution an infinite variance. Painter [11] argued that the increments cannot have an infinite-variance distribution because the
incremental values must be bounded in order to prevent porosity or permeability from taking on nonphysical values, and strictly speaking, this is correct. The truncated PDF have finite variance with power-law decay tails over a large but finite range, resulting in behavior similar to that of the infinite-variance PDF. Therefore, the truncated PDF for lnK increments will be used in the following analysis.

Several approaches have been developed to generate fBm or truncated fLm fields. These methods include the midpoint displacements method, the successive random addition (SRA) method [13], the Fourier filtering method [13], and the modified turning band method [14]. Lu et al. [15] recently compared these methods in detail and concluded that with the SRA algorithm it is easy to understand the geometry and the scaling properties of the stochastic fractals and it is associated with a very attractive computational speed. They developed a simple and efficient three-dimensional simulation code named SRA3DC in Fortran to generate the fBm and the truncated fLm *K* fields. In addition, dispersional and quantile analyses were employed to analyze the increments of the generated fBm and truncated fLm fields to validate the algorithm.

Many researches have been conducted to investigate solute transport in porous media with fBm and / or fLm permeability fields. For example, Painter [6] investigated the statistical moments of travel time for a conservative tracer in the media with an fBm *K* field and a bounded fLm *K* field. Herrick et al. [16] investigated the tritium transport at the highly heterogeneous MADE site; they found that the heavy-tailed *K* fields gave rise to the heavy-tailed velocity fields. Kohlbecker et al. [17] found that a heavy-tailed Levy-stable ln*K* increment field resulted in a heavy-tailed increment field of logarithm velocity (ln v). However, to our knowledge, few researches were about the analysis of plume spatial moment and dispersivity in those media with fLm permeability fields.

In this study, we will use the MCM to investigate solute transport process in the highly heterogeneous aquifer with a truncated fLm field for the lnK increments. The SRA3DC code will be modified and used to generate lnK fields for the MCM analysis. MODFLOW2000 code [18] and MT3DMS code [19] will be used to solve the flow and solute transport problems, respectively. The simulated concentration will be averaged to obtain the ensemble mean concentration for analyzing the plume spatial moments and the longitudinal macro-dispersivity.
## 2. Theoretical Consideration

### 2.1 Lévy-stable Distributions

The Lévy-stable distributions have been studied since the 1930s [20]. Their probability density functions (PDFs) do not have a closed analytic form except for some special cases. Instead, they are expressed by their characteristic function:

$$ \phi(k) = \exp\{-C^{\alpha}|k|^{\alpha}\{1 - i\beta\operatorname{sign}(k)\tan(\alpha\pi/2)\} + iuk\} \quad \alpha \neq 1. \quad (1) $$

where $\alpha$ is the Lévy index with a range of (0,2], $\beta$ is the skewness parameter with a range of [-1,1], $C$ is the width parameter in the interval (0, $\infty$), $\mu$ is the location parameter, $k$ is the Fourier variable, $i^2=-1$, $\tan()$ represents the tangent function and $\operatorname{sign}(k)$ is the sign function of the variable $k$. These characterizing parameters $\alpha$, $C$, $\beta$ and $\mu$ describe the index of stability, the spread, the skewness, and the location of the density, respectively.

When the Lévy-stable distribution is a standard symmetric function with $\beta$ and $\mu$ values equal to zero, Eq. (1) can be simplified as:

$$ \phi(k) = \exp(-C^{\alpha}|k|^{\alpha}). \quad (2) $$

The PDF of symmetric Lévy-stable distribution can be calculated by inverse Fourier transform of Eq. (2)

$$ L_{\alpha}(x) = \frac{1}{\pi} \int_{0}^{\infty} \exp\{-|Ck|^{\alpha} \cos(kx)\} dk. \quad (3) $$

In Eq. (3) the Lévy index $\alpha$ distinguishes the PDFs from one another and characterizes the degree of deviation from the Gaussian distribution which corresponds to the special case $L_{\alpha=2}(x)$. The width parameter $C$ characterizes the discrete extent of variables relative to the mean value, which is similar to the standard deviation of a Gaussian distribution.
## 2.2. Fractional Lévy Motion

The fractional Lévy motion (fLm) is defined as a random and continuous function which has Lévy-stable increments over any lag. The statistical self-affinity or scaling property can be expressed as [21]

$$C_{rh}^{\alpha} = C_{h}^{\alpha} r^{\alpha H} \Rightarrow C_{rh} = r^{H} C_{h} . \quad (4)$$

where $C_{rh}$ and $C_{h}$ are the width parameters for lag distance $rh$ and $h$, respectively. $H$ is the Hurst coefficient, quantifying the degree of interdependence between the increments of fractal process.

When $\alpha=2$, the Lévy-stable distribution becomes a Gaussian distribution, and the fLm model is the same as the fBm model. $\sigma_{rh}$ and $\sigma_{h}$ are the standard deviation of the Gaussian distribution for lag $rh$ and $h$, respectively, corresponding to the width parameters of the Lévy-stable distribution. Substituting $\sigma_{rh}$ and $\sigma_{h}$ into Eq. (4) yields:

$$\sigma_{rh}^{2} = \sigma_{h}^{2} r^{2H} \Rightarrow \sigma_{rh} = r^{H} \sigma_{h} . \quad (5)$$

In general, the Lévy-stable distribution of increments and the scaling property determine the essential features of fLm process. In this paper, we modified Lu's SRA3DC code based on these two features to generate a two-dimensional $\ln K$ field characterized by the truncated fLm distribution. It is necessary to point out that the width parameters $C_1$ and $C_2$ are replaced in the modified code by

$$C_1 = C_0 (1/2)^{H/2} (1 - 2)^{\alpha H/2 - 2} / \alpha . \quad (6)$$

$$C_2 = C_1 (1/2)^{H/2} . \quad (7)$$

where $C_0$ is the given width parameter for random numbers at corner points, $C_1$ and $C_2$ are the width parameters for random numbers at center points and edge middle points, respectively (see Lu et al. [15] for more information).

A truncated fLm realization with $H=0.3$, $C=0.2$ and $\alpha=1.3$ was generated for $\ln K$ field, with which quantile analysis was used to analyze the increments of $\ln K$ to estimate $C$ and $\alpha$. The relationship between the estimated $C$ and the lag is shown in Fig. 1(a). The fitted $H$ value is 0.3259 which is close to the input 0.3 value, and the
approximate straight line on the double logarithmic scale over a wide range is
consistent with fractal model characterized by Eq. (4). In addition, the estimated Lévy
index α equals 1.35 for lag=1 which is slightly larger than the input value of 1.3. On
the other hand, as shown in Fig. 1(b), the Lévy-stable PDF fits the sampled data over
a wider range and can track the power-law decay better into the tails of the
distribution compared with the Gaussian PDF.

<figure><img src="image_2.png" /><figcaption>**Fig. 1 (a).** The relationship between the estimated width parameters and the lag distance; **Fig. 1(b).** Plots of probability density for sample and theoretical Lévy-stable probability density function (PDF) for lag=1</figcaption></figure>

**3. Numerical Simulation**

**3.1 Problem Setup**

A two-dimensional confined aquifer with a domain size 128 m x 128 m was considered with the grid size equaling to 1m (as shown in Fig. 2). Upstream and downstream boundaries were specified as constant head with a head difference of 0.7m in the mean flow direction. No-flow conditions were prescribed at the transverse, upper and lower boundaries. The porosity was considered to be spatially homogeneous with a value of unity for convenience. The geometric mean of lnK was chosen to be 2.5m/d and the increments of lnK were treated as a truncated fLm distribution. The local longitudinal and transverse dispersivities were 1m and 0.1m, respectively. To avoid the boundary effects, a total mass of 1kg was released at (10m, 64m) as a small point source of contamination. And the duration time was finished on the first transport step size.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Schematic illustration of the computational domain</figcaption></figure>

## 3.2 Monte-Carlo Simulation

The MCM was used to simulate solute transport process in the domain. As shown in Table 1, three width parameters, i.e. $C=0.05$ (case 1), 0.1 (case 2) and 0.2 (case 3) were chosen to analyze the impact of width parameters on solute transport process, respectively. Five hundred realizations for each case were carried out in the MCM to satisfy the convergence requirements. The simulation algorithm consists of the following steps:

(1) Using the modified SRA3DC code to generate the random $lnK$ field with its increments following the truncated fLm distribution of given random seed number, mean, $C$, $H$ and $\alpha$ values, respectively. The corresponding $K$ can be obtained with the exponential transformation of the generated $lnK$ random field for each realization.

(2) Solving the flow equation to obtain the velocity field for each realization of the $K$ field. Firstly, flow problem was solved using MODFLOW with the Pre Conditioned Gradient solver (PCG2). Secondly, the cell by cell flow terms were then extracted from the MODFLOW solution and converted from mass flux to velocity.

(3) Solving the transport equation to obtain the concentration field for each realization of the velocity field using MT3DMS with the central finite difference scheme for advection solution. Each realization was stopped when 10% of the initial mass exited the domain.
(4) Repeating steps (1)-(3) for all the realizations and averaging over all realizations to obtain the ensemble mean concentration distribution. The average concentration distribution will be used for moment analysis

<table><caption>**Table 1.** Chatacerizing parameyters considered in this study</caption><thead><tr><th>Case</th><th>α</th><th>C<sub>0</sub></th><th>H</th></tr></thead><tbody><tr><td>1</td><td>1.3</td><td>0.05</td><td>0.3</td></tr><tr><td>2</td><td>1.3</td><td>0.1</td><td>0.3</td></tr><tr><td>3</td><td>1.3</td><td>0.2</td><td>0.3</td></tr></tbody></table>

### 3.3 Moment Analysis

Spatial moments of solute concentration are widely used to analyze solute transport process. The zero-order spatial moment indicating the mass of the system can be expressed as [22]:

$$M_0(t) = \iint_{\Omega} \theta c(x, y, t) dx dy . \qquad (8)$$

where $M_0(t)$ is the mass of the solute in the system at time $t$; $\theta$ is the effective porosity; $c(x, y, t)$ is the concentration at location $(x, y)$ for time $t$, and $\Omega$ represents the interest area.

The first-order moments about the origin are the mass center of solute plume, which can be expressed as:

$$x_c(t) = \frac{1}{M_0} \iint_{\Omega} \theta c(x, y, t) x dx dy . \qquad (9a)$$

$$y_c(t) = \frac{1}{M_0} \iint_{\Omega} \theta c(x, y, t) y dx dy . \qquad (9b)$$

where $x_c(t)$, $y_c(t)$ are the first moments representing the centroid coordinate of solute plume along the longitudinal and transverse directions at time $t$, respectively.

The second-order moments about the plume's center of mass characterizing the macrodispersion can be expressed as:
$$ 
\begin{aligned}
M_{xx}(t) &= \frac{1}{M_0} \iint_{\Omega} \theta c(x, y, t) (x - x_c)^2 dx dy \\
&= \frac{1}{M_0} \iint_{\Omega} \theta c(x, y, t) x^2 dx dy - x_c(t)^2
\end{aligned} 
\quad (10) $$

$$ M_{yy}(t) = \frac{1}{M_0} \iint_{\Omega} \theta c(x, y, t) y^2 dx dy - y_c(t)^2 . \quad (11) $$

where $M_{xx}(t)$, $M_{yy}(t)$ are the second central moments of the concentration plume along the longitudinal and transverse directions at time $t$, respectively.

The longitudinal macrodispersion coefficient $D$ can be expressed as follows:

$$ D = \frac{1}{2} \frac{dM_{xx}}{dt} . \quad (12) $$

Under the assumption of a uniform flow field, the longitudinal macrodispersivity is defined as

$$ \alpha_x(t) = \frac{1}{2} \frac{dM_{xx}}{vdt} . \quad (13) $$

where $\bar{v}$ is the average seepage velocity along the $x$ coordinate, which can be determined by Darcy's law

$$ \bar{v} = \frac{\bar{K} \nabla h}{\theta} . \quad (14) $$

where $\nabla h$ is the hydraulic gradient, $\bar{K}$ is the geometric mean of the permeability, $\theta$ is the porosity which can be treated as a constant.

# 4. Results and Discussions

## 4.1 Permeability Field Properties

Fig. 3 represents realizations of the $\ln K$ fields for case 1-3. It can be found that the spatial distribution patterns of the $\ln K$ fields are similar for different $C$ values. However, the increase in $C$ values leads to the wider disperse range of the $\ln K$ field. Specially, the value of $\ln K$ deviates from 2.5 gradually represented by the enlarged
range of smaller and larger values as C increases. The largest K value is about seven and four orders above its smallest value for C=0.2 and C=0.1, respectively, while the largest K value is one order larger than its smallest value for C=0.05. The results are attributable to the fact that the width parameter C characterizes the dispersion extent of the lnK increments and a larger C value means a larger variation on the lnK increments, which then results in a larger variation in K fields.

<figure><img src="image_4.png" /><figcaption>Fig. 3. Comparison of the lnK field for (1) case 1; (2) case 2; (3) case 3</figcaption></figure>

## 4.2 Velocity Field Properties

Fig. 4 shows the contours of velocity magnitude ($v_x^2+v_y^2)^{1/2}$ in the domain for a realization of case 1-3. It can be found that the increase in C leads to larger variation in the velocity field. The largest value of the velocity magnitude is about 0.12, 0.24 and 0.45 for case 1, 2 and 3, respectively. These results are consistent with the permeability field properties. This is attributed to Darcy's law on condition that the hydraulic gradient and porosity are constant in this study. Consequently, the differences in the permeability field leads to the same trend of the differences in the velocity field.

<figure><img src="image_5.png" /></figure>
<caption>**Fig. 4.** Comparison of the contours of velocity magnitude for (1) case 1; (2) case 2; (3) case 3.</caption>

Combined with the research of Kohlbecker et al. [17] and our present study, it was necessary to make further investigation of the probability distribution of the velocity field. In this study, we mainly examine the impact of lnK with the increments following the truncated fLm distribution on the resulting increments of lnv. For convenience, the longitudinal component velocity vector (i.e. $v_x$) is chosen to be discussed. Fig. 5 is the sampled probability density plot of the increments for lag=1 in lnv$_x$ for different C (i.e. case 2 and 3) values. It indicates that the probability densities of the increments in lnv$_x$ deviate significantly from Gaussian distribution and shift from the body to the tails, which can be approximated well by the Levy-stable distribution. These results are consistent with the results of Kohlbecker et al. that a heavy-tailed distribution of increments in lnK results in a heavy-tailed distribution of increments in lnv. Furthermore, as C increases, the probability densities cover a greater range of values, signifying a larger variation of $v_x$. This is consistent with the variations in the velocity magnitude as shown in Figure 4. In the next section, we will discuss solute transport process obtained from these heavy-tailed velocity fields.

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Comparison of probability density plot of the increments for lag=1 in lnv$_x$ for case 2 and 3</figcaption></figure>

## 4.3 Contour Maps of Solute Concentration

Contour maps of solute concentration are one of the best ways to interpret solute transport process. Fig. 6 shows the comparison of contour maps for case 1-3 at 400 days including a homogeneous case (i.e. C=0) for better analysis. It can be found that solute plume is in Gaussian distribution for C=0 in Fig. 6(a). As C increases, it covers a greater range with a corresponding decrease in peak value, i.e., 7-8 for case 1 in Fig.
6(b), 5-6 for case 2 in Fig. 6(c) and 3-4 for case 3 in Fig. 6(d), respectively. Compared with the homogeneous case, solute plumes for larger C values have a more significant anomalous shape with a sharper leading edge and a wider tailing edge. The concentration gradient decreases with the increase of the C value. All these features reflect more significantly anomalous transport process in a more highly heterogeneous media.

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** The contour maps of the solute concentration for different width parameters *C* at 400 days: (a) *C*=0 (i.e. homogeneous case); (b) case 1; (c) case 2; (d) case 3</figcaption></figure>

**4.4 Spatial Moments**

Fig. 7 is about the temporal variation of the first spatial moments of solute plume for different *C* (i.e. case 2 and 3) values. First of all, the component coordinates of centroids $x_c$ and $y_c$ are basically superposition for different cases through the entire calculation time. Secondly, $x_c$ increases linearly with time while $y_c$ keeps constant for both the cases. These results are consistent with Yan and Wu's study in which the first spatial moments of solute plume are not influenced by the variance of lnK but dominated by the mean of lnK. This is because that the first spatial moments reflecting the mass central position of solute plume are determined by its mean
velocity. As can be seen from Eq. (14), the mean velocity is a linear function of the mean of K with constant hydraulic gradient and porosity. Hence, it was concluded that the first spatial moments are determined by the mean of lnK and independent of the width parameter.

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** The first spatial moments (i.e. $x_c$ and $y_c$) of the solute plume versus time for case 2 and 3</figcaption></figure>

Fig. 8 shows the second spatial longitudinal and transverse moments for different C (i.e. case 2 and 3) values. Fig. 9 shows the longitudinal macrodispersivity $\alpha_x$ obtained using Eq. (13) with the second longitudinal moments from Fig. 8. It's known that for the hypothetical case of homogeneous conditions, with constant velocity and dispersion coefficients and instantaneous point injection, solute plume is Gaussian and the second moments increase linearly with time while the longitudinal macrodispersivity is constant. Our results show that $M_{xx}$ and $M_{yy}$ increase linearly with time on the double logarithmic scale, implying that the second spatial moments can be approximated by a power law function of time, which leads to $\alpha_x$ increase as a power law function of time as shown in Fig. 9. These results may be attributable to the fact that hydraulic properties are random variables characterized by the Lévy-stable distribution, resulting in the spatial distribution of solute concentration deviates from Gaussian distribution and the second spatial moments of solute plume are nonlinear functions of time. This implies that non-Fickian dispersion occurs. Zou et al. [23] have discussed the parabolic relationship between the dispersivity and time, and they have argued that the non-linearly dispersive process is caused by the heterogeneity of the hydraulic properties. Additionally, $M_{xx}$, $M_{yy}$ and $\alpha_x$ increase as C increases, meaning that larger degree of heterogeneity causes larger dispersion of solute plume.
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** The second spatial moments of the solute plume versus time for case 2 and 3: (a) Moments at the longitudinal direction $M_{xx}$; (b) Moments at the transverse direction $M_{yy}$</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Macrodispersivity at the longitudinal direction versus time for case 2 and 3.</figcaption></figure>

**5. Conclusion**

In this study, we used the Monte-Carlo simulation method to investigate the sensitivity of solute transport to the characteristics of the highly heterogeneous and non-stationary random fields. Based on the comparison results as discussed in Section 4, the following conclusions are reached:

Firstly, the heterogeneity of the permeability field is sensitive to the values of the fLm characterizing parameters. A more heterogeneous permeability field is characterized by a larger C value. This is determined by the features of C, which characterizes the discrete extent of variables relative to the mean value and a larger C value leads to a larger variation in K fields.
Secondly, a larger C value lead to a larger variation in the velocity field which are consistent with the permeability field properties. The investigation of the probability distribution shows that the probability density plot of the increments in $lnv_x$ can be approximated well by the Lévy-stable distribution. In addition, as C increases, the probability densities cover a greater range of values, signifying a larger variation of velocity.

Thirdly, contour maps of solute concentration show that larger C leads to more significantly anomalous transport with a sharper leading edge and a wider tailing edge in the plume. The basic reason for this phenomenon is that larger C lead to heavy-tailed K distribution, which directly give rise to heavy-tailed velocity field and then results in more significant anomalous transport.

Fourthly, the first spatial moments of the solute plumes are independent of the fLm characterizing parameters C. The second spatial moments can be approximated by a power law function of time, which leads to the longitudinal macrodispersivities increase as a power law function of time. In addition, both the second spatial moments and the longitudinal macrodispersivities increase as C increases, meaning that a larger degree of heterogeneity causes a larger dispersion of solute plume.

The results summarized above imply that solute transport process is highly dependent on the heterogeneity of the permeability field. The sensitivity analysis underscores the need for careful aquifer characterization for accurate estimation or prediction of subsurface flow and transport. However, our present study only considered two-dimensional case with relatively simple boundary conditions. Further investigation will be conducted for three-dimensional case in highly heterogeneous aquifer with more complex conditions.

## References

1. Dagan, G.: Flow and Transport in Porous Formations. Springer, New York (1989)

2. Hassan, A.E., Cushman, J.H., Delleur, J.W.: Monte Carlo studies of flow and transport in fractal conductivity fields: Comparison with stochastic perturbation theory. Water Resources Research. 33, 2519-2534 (1997)

3. Huang, H., Hassan, A.E., Hu, B.X.: Monte Carlo study of conservative transport in heterogeneous dual-porosity media. Journal of Hydrology. 275, 229-241 (2003)
4. Molz, F.J., Boman, G.K.: A fractal-based stochastic interpolation scheme in subsurface hydrology. Water Resources Research. 29, 3769--377 (1993)

5. Yan, T.T., Wu, J.C.: Impacts of the spatial variation of permeability on the transport fate of solute plume. Advances in Water Science. 17, 30--36 (2006)

6. Painter, S., Mahinthakumar, G.: Prediction uncertainty for tracer migration in random heterogeneities with multifractal character. Advances in Water Resources. 23, 49--57 (1999)

7. Neuman, S.P., Zhang, Y.: A quasi-linear theory of Non-Fickian and Fickian subsurface dispersion.1. Theoretical analysis with application to isotropic media. Water Resources Research. 26, 887--902 (1990)

8. Painter, S.: Random fractal models of heterogeneity: the Lévy-stable approach. *Math Geol*, 27, 813--830 (1995)

9. Liu, H.H., Molz, F.J.: Multifractal analyses of hydraulic conductivity distributions. Water Resources Research. 33, 2483--2488 (1997)

10. Painter, S.: Evidence for non-Gaussian scaling behavior in heterogeneous sedimentary formations. Water Resources Research. 32, 1183--1195 (1996a)

11. Painter, S.: Stochastic interpolation of aquifer properties using fractional Lévy motion. Water Resources Research. 32, 1323--1332 (1996)

12. Meerschaert, M.M., Kozubowski, T.J., Molz, F.J., Lu, S.L.: Fractional Laplace model for hydraulic conductivity. Geophysical Research Letters. 31(8), 1--4 (2008)

13. Peitgen, H.O., Saupe, D.: The science of fractal images. Springer, New York (1988)

14. Yin, Z.: New methods for simulation of fractional Brownian motion. Journal of Computational Physics. 127, 66--72 (1996)

15. Lu, S.L., Molz, F.J., Liu, H.H.: An efficient, three-dimensional, anisotropic, fractional Brownian motion and truncated fractional Lévy motion simulation algorithm based on successive random additions. Computer & Geosciences. 29, 15--25 (2003)

16. Herrick, M.G., Benson, D.A., Meerschaert M.M., McCall K.R.: Hydraulic conductivity, velocity, and the order of the fractional dispersion derivative in a highly heterogeneous system. Water Resources Research. 38, 1227--1239 (2002)

17. Kohlbecker, M.V., Wheatcraft, S.W., Meerschaert, M.M.: Heavy-tailed log permeability distributions imply heavy-tailed log velocity distributions. Water Resources Research. 42, W04411, doi:10.1029/2004WR003815 (2006)

18. McDonald, M.G., Harbaugh, A.W.: A modular three-dimensional finite-difference ground water flow model, USGS Techniques of Water Resources Investigations, Book 6. (1988)

19. Zheng, C.M., Wang, P.P.: MT3DMS: A modular three-dimensional multi-species transport model for simulation of advection, dispersion and chemical reactions of contaminants in
ground water systems: Documentation and user's guide. Contract Report SERDP-99-1, U S
Army Engineer Research and Development Center, Vicksburg, Mississippi. (1999)

20.Lévy, P.: Théorie de l'Addition des Variables Aléatoires. Gauthier-Villars, Paris (1937)

21.Samorodnitsky, G., Taqqu, M.S.: Stable Non-Gaussian Random Processes. Chapman and Hall, New York (1994)

22.Zhang, R.D.: Applied geostatistics in environmental sciences. Science Press, Beijing (2004)

23.Zou, S.M., Xia, J.H., Koussis, A.D.: Analytical solutions to non-Fickian subsurface dispersion in uniform groundwater flow. Journal of Hydrology. 179, 237--258 (1996)