# Vadose Zone Journal

- Soil Science Society of America

## Abstract

It has been recognized that matrix diffusion is an important process for retarding solute transport in fractured rock. Based on analyses of tracer transport data from a number of field tests, we demonstrate for the first time that the effective matrix-diffusion coefficient may be scale dependent and generally increases with test scale. A preliminary theoretical explanation of this scale dependency is also presented, based on the hypothesis that solute travel paths within a fracture network are fractals.

The exchange of solute mass (through molecular diffusion) between fluid in the fractures and fluid in the rock matrix is called matrix diffusion. Owing to the orders-of-magnitude slower flow velocity in the matrix compared with fractures, matrix diffusion can significantly retard solute transport in fractured rock, and therefore is an important process for a variety of problems, including remediation of subsurface contamination and geological disposal of nuclear waste (e.g., Neretnieks, 2002; Jardine et al., 1999; Bodvarsson et al., 2000; Guimera and Carrera, 2000).

The effective matrix diffusion coefficient (molecular diffusion coefficient in free water multiplied by matrix tortuosity) is an important parameter for describing matrix diffusion and in many cases largely determines overall solute transport behavior. For example, Fig. 1 shows simulated breakthrough curves for a radionuclide moving from a repository to the water table within the unsaturated zone of Yucca Mountain, Nevada, the proposed site for a high level nuclear waste repository. These simulated curves were generated with the site-scale model of Yucca Mountain (Wu et al., 2002; Moridis et al., 2003). The distance between the repository horizon and the water table is about 300 m. As shown in Fig. 1, the performance of the proposed repository (or travel time of a radionuclide within the unsaturated zone) is largely determined by the effective matrix-diffusion coefficient. A larger coefficient corresponds to a better performance.

Matrix diffusion coefficient values measured from small rock samples in the laboratory are generally used for modeling field-scale solute transport in fractured rock (Boving and Grathwohl, 2001; Moridis et al., 2003). 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 (Shapiro, 2001; Neretnieks, 2002; Liu et al. 2003a, 2003b). These authors have also offered several different explanations for the need to increase coefficient values at test sites under consideration. However, the fundamental question of whether a systematic relation exists between effective matrix diffusion coefficient and test scale remains unanswered.

The major objective of this note is to demonstrate, based on data published in the literature, that the effective matrix-diffusion coefficient may be scale dependent and generally increases with scale. A preliminary explanation, based on a hypothesis of fractal paths for solute transport in fractures, is presented for this possible scale-dependent behavior.

## Data from Literature

Effective matrix diffusion–coefficient values have been estimated from a number of field test sites characterized by different rock types. To compile these values (corresponding to different tracers) as a function of test scale, we define a matrix diffusion–coefficient ratio, RD, as an effective coefficient value (estimated from field data) divided by a local value. It is an indicator of scale dependency that is expected to exist when RD is always larger than one at field scale and is a function of the test scale.

The local matrix diffusion–coefficient values, *D*_{e}, refer to laboratory measurements for small rock samples or estimations from the following relationship (when laboratory measurements are not available): 1where *D*_{0} is the molecular diffusion coefficient in free water, and τ is the tortuosity factor, determined here with Archie's Law (Boving and Grathwohl, 2001): 2Here, φ is the matrix porosity, and *m* is an empirical parameter. Boving and Grathwohl (2001) compiled *m* values for different types of rocks and found that *m* is generally >2 in materials of low porosity (≤0.2). By comparing measured local-scale matrix diffusion coefficients and the corresponding molecular diffusion–coefficient values in free water, Moridis and Hu (2000) concluded that *m* = 2 is a good approximation for the rock matrix in the unsaturated zone of Yucca Mountain. The laboratory data cited by Becker and Shapiro (2000) suggest that *m* = 2.93 (obtained using measured φ = 1.5%) is valid for the crystalline rock matrix under consideration. Therefore, it is very likely that *m* = 2 corresponds to the upper limit for the tortuosity factor for a rock matrix with porosity <0.2. To avoid potential exaggeration of scale effects (or an artificial increase in estimated RD values), we use *m* = 2 here. Note that some researchers (Maloszewski and Zuber, 1993; Jardine et al., 1999) also used effective diffusion coefficients estimated from Eq. [1] to match field observations. In this case, the corresponding RD values are simply determined as the ratio of their selected values for tortuosity factor to the corresponding rock matrix porosity (≤0.2).

Figure 2 shows the relationship between effective matrix diffusion coefficients, determined from a number of sites by different research groups (Maloszewski and Zuber, 1993; Jardine et al., 1999; Becker and Shapiro, 2000; Callahan et al., 2000; Shapiro, 2001; Neretnieks, 2002; Liu et al., 2003a, 2003b) and the corresponding test scales. Jardine et al. (1999) performed and analyzed a carefully designed tracer test in fractured shale bedrock. They were able to obtain breakthrough curves from wells with different distances to the source. However, the same parameters were used to match all the breakthrough curves for a given tracer. As indicated in their Fig. 9, the breakthrough curves from wells about 6 m from the source are reasonably matched, but the curve from the well 17 m from the source is not. A larger effective matrix diffusion–coefficient value seems to be needed to match the data from the latter well, an indication of a possible increase in the coefficient with scale. Therefore, in Fig. 2, the fitted effective matrix diffusion is considered to correspond to a test scale of 6 m only. Note that Jardine et al. (1999) indicated their intentions were not to rigorously model all of the processes contributing to solute transport in the system, but rather to test the importance of matrix diffusion.

Neretnieks (2002) reported matches to tracer test data collected from the Äspö site with a test scale of 5 m and found a need for a factor 30 times larger for the fracture–matrix interface area (or effective matrix-diffusion coefficient) than expected. Note that the increase in fracture–matrix interface area is equivalent to the increase in effective diffusion coefficient (for a given interface area in a model). Interestingly, he also indicated that nine other research groups had also independently evaluated the tracer test data from the site using different modeling approaches. Nearly all the groups found the need for a factor 30 to 50 times larger effective fracture–matrix interface area (or effective matrix-diffusion coefficient) than expected. In Fig. 2, a representative RD value of 40 is used for the Äspö test site.

Liu et al. (2003a)(2003b) presented model analyses of two different field test data, collected in the unsaturated zone of Yucca Mountain by scientists from the U.S. Geologic Survey and Lawrence Berkeley National Laboratory. Unlike studies reported by other researchers mentioned in this note, Liu et al. (2003a)(2003b) matched both the flow field (characterized by water travel time and/or seepage into subsurface openings) and tracer breakthrough curves. They reported that increased fracture–matrix interface areas (or effective matrix diffusion coefficients) were needed for both tests. Note that the data of Callahan et al. (2000) were also collected for the rock matrix in the unsaturated zone of Yucca Mountain.

Becker and Shapiro (2000) and Shapiro (2001) reported analyses of tracer test data from fractured crystalline rock at Mirror Lake site. Becker and Shapiro (2000) showed that laboratory measurement of the effective diffusion coefficient should be replaced by the coefficient in free water (corresponding to RD = 3333) to match the bromide data in their Test C with a test scale of about 36 m. However, they were not able to match all the breakthrough curves for different tracers and argued that advective transport processes contribute to this discrepancy. An alternative explanation may be that a simple model used by those authors cannot capture all the importance processes, such as effects of subsurface heterogeneity, even when matrix diffusion is a dominant process. Nevertheless, the value of RD = 3333 is included in Fig. 2. Shapiro (2001) found that three to five orders of magnitude greater than the estimates of the matrix-diffusion coefficient from laboratory experiments were needed to match the tracer data observed at a kilometer scale. His analysis probably provides the first estimate for kilometer-scale effective diffusion coefficient. A representative value of Shapiro (2001) for RD (1 × 10^{4}) is used in Fig. 2.

Figure 2 also includes analyses of results of Maloszewski and Zuber (1993) for three different sites. Note that a recent analysis result of a tracer test by McKenna et al. (2001) is not included in Fig. 2 because the analysis does not provide an independent estimate of the effective diffusion coefficient, but rather a lump parameter that includes effective diffusion coefficient, fracture spacing, and other properties.

## A Fractal-Based Explanation

Although some uncertainties exist, the data shown in Fig. 2 seem to strongly suggest that the effective matrix diffusion coefficient, like permeability and dispersivity, increases with test scale. A number of researchers have attempted to explain why the effective coefficient determined from field data is larger than the corresponding laboratory value (Shapiro, 2001; Neretnieks, 2002; Liu et al., 2003a, 2003b).

Shapiro (2001) suggested that kilometer-scale “effective matrix diffusion” is not a diffusive process, but actually an advective process between high and low permeability zones, resulting in a significantly large “effective diffusion coefficient.” While this may be a plausible explanation, further confirmation is still needed. For example, Liu et al. (2003a)(2003b) used a dual-permeability model involving both fast flow in fractures and slow flow in the matrix (as well as the advective transport between the two) and still found the need to use increased effective diffusion coefficients for matching the tracer test data. Neretnieks (2002) argued that the existence of fracture in-filling creates relatively large areas for solute to diffuse into rock matrix, which, together with the process of diffusion into stagnant water, contributes to the need for increasing the effective diffusion coefficient to match the data. Liu et al. (2002)(2003a, 2003b) and Wu et al. (2003) indicated that the existence of many small-scale fractures (which considerably increase the fracture–matrix interface area but are not considered in numerical models) may be the major reason for the relatively large effective diffusion coefficient calculated from field data. While these suggested mechanisms seem to be reasonable for field-scale solute transport in fractured rock, they cannot be directly used to explain why the effective diffusion coefficient increases with test scale.

In this study, we propose a fractal-based explanation for the scale-dependent behavior of the effective diffusion coefficient. The fractal concept has been found to be useful for describing both subsurface heterogeneity and many flow and transport processes (e.g., Wheatcraft and Tyler, 1988; Molz and Boman, 1993). In commonly used numerical and analytical models of solute transport including matrix diffusion (e.g., Sudicky and Frind, 1982; Wu et al., 2003), an actual fracture network is generally conceptualized as parallel vertical or horizontal fractures. A fracture wall is approximated as a flat wall. In this case, solute particle travel paths within fractures are generally straight lines. However, the actual solute particle travel path is much more intricate and tortuous for the following reasons. First, fracture walls are not flat but rough. The rough surface generates a much larger fracture–matrix interface area than a flat fracture wall, and the fracture roughness is characterized by fractals (National Research Council, 1996). Second, fractures exist at different scales, with small-scale fractures generally excluded from modeling studies (Wu et al., 2003). However, these small-scale fractures can make flow and transport paths much more tortuous than straight lines, as demonstrated by Liu et al. (2002). Considering that both fracture roughness and fracture-network geometry can be characterized by fractals (e.g., Barton and Larsen, 1985), it is reasonable to hypothesize that a solute travel path within a fracture network is fractal, rather than a straight line, as is assumed in many numerical or analytic models.

The length of the fractal solute travel path (*L*) between tracer release and monitoring points depends on the spatial scale (or length of a ruler, δ) used to measure it. We denote the straight-line distance between the release and monitoring points as *L**. Considering that matrix diffusive flux is proportional to the fracture–matrix interface area (or the length of solute travel path), and assuming the local matrix-diffusion coefficient to be precise, we can approximate RD as a ratio of actual fracture–matrix interface area to the area used in numerical or analytical models: 3

Following the procedure of Feder (1988) to determine the length of the coast of Norway, *L* can be determined by 4where *N* is the number of rulers with length δ needed to measure the length of fractal solute travel path between tracer release and monitoring points, and *D* > 1 is the fractal dimension of a solute travel path. Assuming that a solute travel path within a small interval δ can be approximated as a straight line, we obtain the following relationship based on Eq. [3] and [4]: 5The above equation indicates that RD is a power function of test scale *L**. Because *D* > 1, RD increases with *L**, consistent with results showed in Fig. 2. The power value, *D* − 1, in Eq. [5] is expected to be site specific. Fitting Eq. [5] to data points (Fig. 2) for the Yucca Mountain site (Callahan et al., 2000; Liu et al., 2003a, 2003b) and Mirror Lake site (Becker and Shapiro, 2000; Shapiro, 2001) results in *D* = 1.7 and 1.3, respectively.

## Conclusions

While the scale dependency of permeability and dispersivity has been a popular research topic for many years in the subsurface hydrology community (e.g., Gelhar, 1993; Neuman, 1994), we, for the first time, demonstrate in this note that the effective matrix-diffusion coefficient may also be scale-dependent and increases with test scale. This finding has many important implications for problems involving matrix diffusion. As clearly demonstrated in Fig. 1, the simulated radionuclide travel time within the unsaturated zone of Yucca Mountain may be significantly underestimated when this scale-dependent behavior is not considered.

We also presented a fractal-based explanation for the possible scale-dependent behavior of the effective matrix-diffusion coefficient. However, uncertainties exist in the estimated effective diffusion coefficients given in Fig. 2 because these coefficients have been obtained from inverse modeling, which cannot give unique parameter values. Therefore, more carefully designed field tests and numerical experiments using realistic fracture networks are still needed to confirm this scale-dependent behavior and to develop more rigorous theoretical explanations. Furthermore, non-Gaussian solute transport behavior, as discussed by a number of authors (e.g., Berkowitz and Scher, 1998; Painter et al., 2002), has not been used in deriving the effective matrix diffusion coefficient values shown in Fig. 2. How the consideration of this non-Gaussian behavior affects the observed scale-dependency of matrix diffusion coefficient also needs to be investigated in the future.

## Acknowledgments

We are indebted to Q. Zhou and D. Hawkes at Lawrence Berkeley National Laboratory for their careful review of a preliminary version of this manuscript. We also thank two anonymous reviewers for their constructive comments. This work was supported by the Director, Office of Civilian Radioactive Waste Management, U.S. Department of Energy, through Memorandum Purchase Order EA9013MC5X between Bechtel SAIC Company, LLC, and the Ernest Orlando Lawrence Berkeley National Laboratory (Berkeley Lab). The support is provided to Berkeley Lab through the U.S. Department of Energy Contract no. DE-AC03-76SF00098.

- Received July 10, 2003.