Published online 23 August 2007
Published in Vadose Zone J 6:679-683 (2007)
DOI: 10.2136/vzj2006.0140
© 2007 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
TECHNICAL NOTES
Scale Dependence of the Effective Matrix Diffusion Coefficient: Some Analytical Results
H. H. Liua,*,
Y. Q. Zhanga and
F. J. Molzb
a Earth Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA
b Dep. of Environmental Engineering and Science, Clemson Univ., Clemson, SC
* Corresponding author (hhliu{at}lbl.gov).
All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.
ABSTRACT
Matrix diffusion is an important process affecting solute transport in fractured rock, and the matrix diffusion coefficient is a key parameter for describing this process. Previous studies have indicated that effective matrix diffusion coefficient values obtained from a number of field tracer tests are enhanced in comparison with local values and may increase with test scale. In this communication, we develop analytical expressions for the effective matrix diffusion coefficient for two simple fracture–matrix systems and demonstrate that heterogeneities in the rock matrix at different scales may contribute to the scale dependence of the effective matrix diffusion coefficient.
The matrix diffusion coefficient is an important parameter relating to transport in fractured rock, one that in many cases largely determines overall solute transport behavior. Recently, several research groups have independently found that effective matrix diffusion coefficients much larger than laboratory measurements are needed to match field-scale tracer-test data (e.g., Neretnieks, 2002; Becker and Shapiro, 2000; Shapiro, 2001; Liu et al., 2004a). These observations imply that the effective matrix diffusion coefficient may be scale dependent and increase with test scale.
In addition to the observed enhancement, Liu et al. (2004b) and Zhou et al. (2007) examined many field-scale tracer-testing results available in the literature and found that scale dependence of the matrix diffusion coefficient was highly likely. This potential scale dependence has important implications for large-scale solute transport in fractured rock. Liu et al. (2007) used numerical experiments to show that a combination of local flow loops (formed by locally connected small-scale fractures) and the associated matrix diffusion process, together with scaling properties in flow-path geometry, can cause the effective (or apparent) matrix diffusion coefficient to increase with travel distance, similar to the well-known scale dependence of dispersivity.
Given the complex nature of fracture–matrix systems, it is likely that more than flow-path geometry alone contributes to scale dependence. The major objective of this note is to demonstrate that rock matrix heterogeneities may also contribute to the scale dependence of the effective matrix diffusion coefficient. In particular, we develop analytical expressions for the effective matrix diffusion coefficient for two idealized fracture–matrix systems because analytical results can generally provide useful physical insights and are transparent.
Effective Matrix Diffusion Coefficient for a Single Fracture
Analytical solutions for solute transport (with matrix diffusion) in a homogeneous system involving a single fracture are commonly used to analyze and interpret field observations (Neretnieks, 2002). Here we determine the effective matrix diffusion coefficient for a single fracture involving heterogeneous rock matrix properties (Fig. 1
).
View larger version (11K):
[in this window]
[in a new window]
|
FIG. 1. Schematic of a single fracture system. The parameter b refers to the half aperture of the fracture.
|
|
In our derivation, several assumptions (or simplifications) are made. Following previous studies (Tang et al., 1981; Neretnieks, 2002), we assume fracture aperture and water velocity to be constants and advection between fracture and matrix to be negligible. The neglect of the advection is justified because in general, rock matrix permeability is significantly smaller than the corresponding fracture permeability. However, this does not mean that rock matrix permeability and its spatial variability do not exist for the given problem under consideration. To explore the effects of rock matrix heterogeneity, we assume rock matrix to be characterized by a layered system. This kind of treatment has been often used as a first step to deal with interaction between flow and transport processes and subsurface heterogeneity (Güven et al., 1984; Yeh and Harvey, 1990). While more studies are needed to address heterogeneity effects in more general cases, the layered-system assumption can be justified roughly by field test conditions in a fractured rock. For most of these tests, the tracer penetration depth due to matrix diffusion is significantly smaller than the test scale (or tracer transport length along the fracture). As a result, the degree of matrix heterogeneity along the fracture direction tends to be much larger than that normal to the fracture (within the range of penetration depth). We also assume a perfect correlation between rock permeability and porosity. This allows us to develop analytical results, while we expect that consideration of a more realistic correlation would not change our conclusions qualitatively regarding the effects of rock matrix heterogeneity.
Our derivation is based on recent theories for determining solute transport in a single fracture with a wide range of retention processes (including matrix diffusion) and spatially variable flow and transport properties (Cvetkovic et al., 2004; Painter and Cvetkovic, 2005). According to these theories, the impulse–response function in the time domain for such a single fracture system with unlimited fracture spacing, which also may be viewed as the probability density distribution for a unit pulse input of conservative solute, is given by (Painter and Cvetkovic, 2005)
 | [1] |
where H is the Heaviside function, and t is time. The residence time 
is defined by
 | [2] |
where l is the distance between the inlet and the location where a breakthrough curve is observed, and V is the water flow velocity along a fracture. Thus, l can be considered to represent the "scale" of the system. The parameter B is defined as
 | [3] |
where 
, D, and b are the matrix porosity, local matrix diffusion coefficient (molecular diffusion coefficient multiplied by tortuosity factor), and local half aperture, respectively. Note that the form of Eq. [3] is slightly different from the definition of parameter B given in Painter and Cvetkovic (2005) because of the difference between definitions of diffusion coefficients. The diffusion coefficient in Painter and Cvetkovic (2005) is the same as D multiplied by . The results based on Eq. [1] and [3] are consistent with the analytical solution for solute transport in a single fracture given by Neretnieks (2002). On the basis of a sizable number of laboratory measurements, Boving and Grathwohl (2001) showed that the tortuosity factor could be approximated by the so-called Archie's law (a power function of the matrix porosity) for different types of rocks. This allows D to be further expressed as
 | [4] |
where D0 is the molecular diffusion coefficient in free water, and the empirical parameter n is generally >1 in materials of low porosity (<0.2) (Boving and Grathwohl, 2001). Similarly, for a heterogeneous system, a representative local-scale matrix diffusion coefficient is given by
 | [5] |
where m is the mean matrix porosity (corresponding to the geometric mean for a log-normal porosity distribution of local porosities). Note that D can be spatially variable as a result of matrix porosity heterogeneity, and Dm is a constant for a given problem.
For simplicity, we consider a fracture–matrix system (Fig. 1) with uniform matrix porosity values in the z direction but spatially variable values in the x direction (a layered system). Based on Eq. [3], and keeping in mind that effective parameters are obtained by replacing the real, heterogeneous system with an effective, homogeneous system, we have
 | [6] |
where D* is the effective matrix diffusion coefficient for the test scale l. Combining Eq. [3] and [6], and considering bV to be constant, yields
 | [7] |
Consider to be a random variable that has a log-normal distribution with a standard deviation of 
l (for the log of the matrix porosity) within the interval between x = 0 and l. (The reasoning for choosing the log-normal distribution will be discussed later.) Using the well-known expression for the arithmetic mean and moments of a log-normally distributed property and noting that 1 + n/2 also has a log-normal distribution (because is log-normally distributed), we have
 | [8] |
where 
refers to a mean for a spatially variable parameter between x = 0 and l, and ln(1+n/2) refers to the standard deviation of ln(1 + n/2) between x = 0 and l. Based on Eq. [8] and using the following properties:
 | [9] |
 | [10] |
 | [11] |
Eq. [7] can be further expressed as
 | [12] |
If the log of rock matrix porosity is characterized by a stationary stochastic process, l is not dependent on the scale l. However, a number of studies indicate that subsurface heterogeneity follows fractal-like long-range correlations that will give a scale-dependent l (e.g., Molz et al., 1997). Note that the above derivation is based on the ergodicity approximation that allows for estimating meaningful heterogeneity statistics from a single realization.
While studies on porosity variability are relatively limited (Hassan et al., 1998), there are many studies in the literature on the permeability (k) variation and scaling of ln(k). Therefore, it is of interest to relate the potential scale dependence of D* to that for permeability. For the given rock matrix shown in Fig. 1, the effective matrix permeability (for water flow in the z direction) is the same as the arithmetic mean of permeability values between x = 0 and l. Like many other researchers, we consider permeability as following log-normal distribution with a standard deviation F for ln(k). In this case, using the well-known expression for the arithmetic mean and moments of a log-normally distributed property, the effective permeability (ke) is given as
 | [13] |
where kg is the geometric mean of local permeability values.
Different relations between porous-medium permeability and porosity have been proposed in the literature. Such a relation was recently provided by Costa (2006):
 | [14] |
In his derivation of Eq. [14], Costa (2006) considered the effects of pore size on rock permeability, in addition to other factors. Because a rock matrix generally has small porosity values, the above equation can be further simplified to
 | [15] |
The above equation is also supported by study results of Hassan (2001). He analyzed measured data for porosity and conductivity (proportional to permeability) collected from two test sites—Amchitka Island, AK, and the Central Nevada Testing Area—and found that the conductivity could be approximated by a power function of porosity. The assumption of log-normal distribution for matrix porosity can be justified from Eq. [15] as a direct result of log-normal distribution for permeability. The same assumption was also used by Hassan et al. (1998) to explore the effects of porosity heterogeneity on solute transport in porous media.
Combining Eq. [12], [13], and [15] yields
 | [16] |
where m and n are exponents characterizing the porosity–permeability relationship and tortuosity for matrix diffusion, respectively. Assuming a typical value of m = 3 (Costa, 2006) and a value of n = 1 corresponding to the lower limit of n (Boving and Grathwohl, 2001; Liu et al., 2004b), the exponent in the Eq. [16] is equal to 0.5.
Equations [12] and [16] establish the intrinsic relations among effective matrix diffusion, effective matrix permeability, and porosity variability. It has been widely recognized that effective permeability is scale dependent and generally increases with test scale; Eq. [16] indicates that the effective matrix diffusion coefficient should follow the same trend as the permeability.
Effective Matrix Diffusion Coefficient for Multiple Flow Channels
Whereas the focus of the previous section was on solute transport in a single fracture, flow and transport processes in fractured rock are in fact characterized by different flow channels within a fracture network (e.g., Neretnieks, 2002). Different flow channels may have different flow and transport properties. In this study, we use a simplified conceptual flow model to investigate the effects of interchannel heterogeneity of diffusive properties. Specifically, we consider a simplified multichannel system in which each flow channel has uniform properties and does not mix with any other channels, except at influent and effluent points. Other researchers have used similar conceptual models to analyze flow and transport processes in fractured rock (e.g., Neretnieks, 2002). These channels have the same length, width, fracture aperture, and water velocity but different matrix diffusion coefficients. This treatment allows us to exclude the effects of dispersion processes and focus on the effects of variability in diffusive properties among different flow channels. Following other researchers (Neretnieks, 2002; Painter and Cvetkovic, 2005), we assume unlimited fracture spacing for calculating solute transport along each flow channel. This is justified by the fact that for most field-scale tests, diffusion–penetration depths (characterized by 
) are much smaller than fracture spacing. When the above condition is not satisfied, our derivation to be presented here will cease to be valid.
We define parameter a as 
and assume a to follow a normal distribution.When t > = x/V and for a pulse input (with mass M0), the solute concentration (c) for a channel is given by (Tang et al., 1981)
 | [17] |
For t < , c(x, t) = 0 (Tang et al., 1981). In the original development of Tang et al. (1981), the analytical solution was presented for step-pulse tracer release at the fracture inlet. However, it can be shown that Eq. [17] is identical to Tang et al. (1981) when a pulse input is considered. We assume the probability density function for a to be
 | [18] |
where is the standard deviation of a, and 
is the mean of a. So, the average concentration (for t > ) at the effluent point is (Zhang et al., 2006)
 | [19] |
The normal distribution for a is assumed here simply because it allows the derivation of an analytical solution for the average concentration, while the actual distribution for a given field-scale problem is not well understood at this stage. However, this analytical result is valid only for a small degree of heterogeneity, such that an insignificant portion of negative a values is included in Eq. [18]. To do so, we may impose the following limit to Eq. [18] and [19]:
 | [20] |
Equation [19] can be further written as
 | [21] |
where
 | [22] |
 | [23] |
The analogy between Eq. [17] and [21] allows us to determine effective transport parameter values from observed breakthrough curves resulting from a combination of multiple flow channels. Specifically, a comparison of Eq. [17] and [21] indicates that Eq. [17] can be used to fit the average breakthrough curve (Eq. [21]) perfectly when T and A are treated as effective values for and a, respectively. (Note that Eq. [22] and [23] are valid only when the denominator in [23] is positive.) In this case, the ratio of the effective matrix diffusion coefficient (Fd) to a representative local-scale coefficient is
 | [24] |
For a given velocity, the residence time is proportional to the distance (or test scale). For typical values of V = 1 m d–1 and = 0.001 s–0.5 (Zhou et al., 2007) and for ' = 1/3, Fig. 2
shows Fd as a function of test scale. Again, this analytical result indicates the scale dependence of the effective matrix diffusion coefficient. Note that for a given test scale, an estimated effective matrix diffusion coefficient depends on velocity V. A large velocity would give a smaller residence time and therefore a smaller effective matrix diffusion coefficient (see Eq. [23]).
The effective residence time T has a complex relation with the actual residence time . It initially increases and then decreases with . In their numerical studies, Zhang et al. (2006) could not detect the scale dependence of the effective matrix diffusion coefficient for a fracture system with interchannel heterogeneity, because they assumed T = . They also noted that their matches with the average breakthrough curves were not satisfactory and limited their simulations to those with relatively small heterogeneities.
Summary and Conclusions
By developing analytical results for two fracture–matrix systems, we have demonstrated here that rock matrix heterogeneity is likely to contribute to the scale dependence of the effective matrix diffusion coefficient, with the effective coefficient increasing with scale. These analytical results can provide insightful relations between the effective matrix diffusion coefficient and other parameters such as porosity and hydraulic conductivity.
We derived relationships among the effective matrix diffusion coefficient, rock matrix porosity variability, and effective rock matrix permeability for a single fracture system with spatially variable local matrix diffusion coefficient. These relationships clearly indicate that the effective matrix diffusion coefficient and rock matrix permeability follow a similar trend in terms of their relationship to test scale, with permeability recognized in the literature as increasing with test scale. We also developed an analytical expression for the effective matrix diffusion coefficient in a fracture–matrix system involving multiple flow channels, each of which has its own matrix diffusion coefficient. For a given water-flow velocity, the expression again suggests that the effective matrix diffusion coefficient keeps increasing with test scale, as a result of the variability in diffusive properties among different flow channels. This solution also indicates that for a given test scale, a higher water velocity yields a smaller effective matrix diffusion coefficient. However, we emphasize that these analytical results were developed for the two idealized systems to demonstrate the potential relationship between subsurface heterogeneity at different scales and the observed scale dependence of the effective matrix diffusion coefficient. More studies are needed to determine analytical relations regarding heterogeneity effects in more general field-scale cases.
In natural systems, the scale dependence of the effective matrix diffusion coefficient likely results from a combination of a number of different mechanisms. So far, we have identified two important mechanisms: the complexity of fracture-network geometry, which has been ignored in our current modeling practice (Liu et al., 2007), and the different-scale subsurface heterogeneity discussed in this note. Both of these neglected mechanisms will be important in field-scale simulations of radionuclide and other solute transport in fracture–matrix systems.
ACKNOWLEDGMENTS
We are indebted to Jim Houseworth and Daniel Hawkes at Lawrence Berkeley National Laboratory (LBNL) for their careful review of a preliminary version of this manuscript. We also appreciate the critical and constructive review comments from Dr. Aimo Hautojarvi from Posiva, Finland, and five anonymous reviewers. This report was prepared by LBNL pursuant to a contract funded by the USDOE, Office of Civilian Radioactive Waste Management (OCRWM), Office of Science and Technology and International (OST&I), and neither LBNL nor any of its contractors or subcontractors nor the USDOE/OCRWM/OST&I nor any person acting on behalf of either makes any warranty or representation, express or implied, with respect to the accuracy, completeness, or usefulness of the information contained in this report; or that the use of any information, apparatus, method, or process disclosed in this report may not infringe privately owned rights; or assumes any liabilities with respect to the use of, or for damages resulting from the use of, any information, apparatus, method, or process disclosed in this report. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by USDOE/OCRWM/OST&I. The views, opinions, findings, and conclusions or recommendations of authors expressed herein do not necessarily state or reflect those of the USDOE/OCRWM/OST&I.
REFERENCES
- Becker, M.W., and A.M. Shapiro. 2000. Tracer transport in fractured crystalline rock: Evidence of non-diffusive breakthrough tailing. Water Resour. Res. 36(7):1677–1686.[CrossRef]
- Boving, T.B., and P. Grathwohl. 2001. Tracer diffusion coefficients in sedimentary rock: Correlation to porosity and hydraulic conductivity. J. Contam. Hydrol. 53:85–100.[CrossRef][Web of Science][Medline]
- Costa, A. 2006. Permeability–porosity relationship: A reexamination of the Kozeny–Carman equation based on a fractal pore-space geometry assumption. Geophys. Res. Lett. 33:L02318, doi:10.1029/2005GL025134.[CrossRef]
- Cvetkovic, V., S. Painter, N. Outters, and J.-O. Selroos. 2004. Stochastic simulation of radionuclide migration in discretely fractured rock near Äspö Hard Rock Laboratory. Water Resour. Res. 40:W02404, doi:10.1029/2003WR002655.[CrossRef]
- Güven, O., F.J. Molz, and J.G. Melville. 1984. An analysis of dispersion in a stratified aquifer. Water Resour. Res. 20(10):1337–1354.
- Hassan, A.E. 2001. Water flow and solute mass flux in heterogeneous porous formations with spatially random porosity. J. Hydrol. 242:1–25.[CrossRef]
- Hassan, A.E., J.H. Cushman, and J.W. Delleur. 1998. Significance of porosity variability to transport in heterogeneous porous media. Water Resour. Res. 34(9):2249–2259.[CrossRef]
- Liu, H.H., R. Salve, J.S. Wang, G.S. Bodvarsson, and D. Hudson. 2004a. Field investigation into unsaturated flow and transport in a fault: Model analyses. J. Contam. Hydrol. 74(1–4):39–59.[CrossRef][Web of Science][Medline]
- Liu, H.H., G.S. Bodvarsson, and G. Zhang. 2004b. Scale dependency of the effective matrix diffusion coefficient. Vadose Zone J. 3:312–315.[Abstract/Free Full Text]
- Liu, H.H., Y.Q. Zhang, Q. Zhou, and F.J. Molz. 2007. An interpretation of potential scale dependence of the effective matrix diffusion coefficient. J. Contam. Hydrol. 90:41–57.[CrossRef][Web of Science][Medline]
- Molz, F.J., H.H. Liu, and J. Szulga. 1997. Fractional Brownian motion (fBm) and fractional Gaussian noise (fGn) in subsurface hydrology. Water Resour. Res. 33(10):2273–2286.[CrossRef]
- Neretnieks, I. 2002. A stochastic multi-channel model for solute transport: Analysis of tracer tests in fractured rock. J. Contam. Hydrol. 55:175–211.[CrossRef][Web of Science][Medline]
- Painter, S., and V. Cvetkovic. 2005. Upscaling discrete fracture network simulations: An alternative to continuum transport models. Water Resour. Res. 41:W02002, doi:10.1029/2004WR003682.[CrossRef]
- Shapiro, A.M. 2001. Effective matrix diffusion in kilometer-scale transport in fractured crystalline rock. Water Resour. Res. 37(3):507–522.[CrossRef]
- Tang, D.H., E.O. Frind, and E.A. Sudicky. 1981. Contaminant transport in a fractured porous media: Analytical solution for a single fracture. Water Resour. Res. 17(3):555–564.
- Yeh, T.-C.J., and D.J. Harvey. 1990. Effective unsaturated hydraulic conductivity of layered sands. Water Resour. Res. 26(6):1271–1279.[CrossRef]
- Zhang, Y.Q., H.H. Liu, Q. Zhou, and S. Finsterle. 2006. Effects of dual-scale heterogeneity on effective matrix diffusion coefficient in fractured rocks. Water Resour. Res. 42:W04405, doi:10.1029/2005WR004513.[CrossRef]
- Zhou, Q., H.H. Liu, F.J. Molz, Y. Zhang, and G.S. Bodvarsson. 2007. Effective matrix diffusion coefficient for fractured rock: Results from a literature survey. J. Contam. Hydrol. (in press).
