Numerical Analysis of a Short-Term Tracer Experiment in Fractured Sandstone

1 Department of Earth and Environmental Sciences, National Chung Cheng University, Chiayi, Taiwan, ROC * Corresponding author address: Prof. Tai-Sheng Liou, Department of Earth and Environmental Sciences, National Chung Cheng University, Chiayi, Taiwan, ROC; E-mail: seiliou@eq.ccu.edu.tw doi: 10.3319/TAO.2007.18.5.1029(Hy) A short-term, pulse injection tracer experiment conducted in fractured quartzitic sandstone at Kukuan, Taiwan was analyzed. Tracer transport at the test site was dominated by advection but a specific attenuation mechanism leading to breakthrough curve (BTC) tailing also seemed to exist. Matrix diffusion was hypothesized as the transport mechanism that results in the tailing. This hypothesis was proved by comparing the field BTC with numerical simulation results obtained by the general-purpose flow/transport simulator, TOUGH2, based on a single-fracture conceptual model. Due to the lack of accuracy of estimating the interporosity flux by the conventional double porosity model (DPM), TOUGH2 was incorporated with the multiple interacting continua (MINC) scheme to simulate the transient characteristics of the interporosity flux. In MINC, rock matrix is discretized as a series of continua according to the perpendicular distance from the fracture that adjoins the matrix. The closer the rock matrix is to the fracture, the finer the rock matrix is discretized. This concept is fundamentally different from DPM in that rock matrix is no longer treated as a single continuum. Simulation results by TOUGH2-MINC have successfully reproduced the observed BTC tailing even under the dominating advection effect. Sensitivity studies showed that TOUGH2-MINC is sensitive to parameters including fracture aperture (2b), matrix porosity (nm) and effective molecular diffusion coefficient in matrix (Dm). If 2b, nm , Dm , are respectively 200 μm, 2%, 10 m s , and if hydrodynamic dispersion coefficient Terr. Atmos. Ocean. Sci., Vol. 18, No. 5, December 2007 1030 (D) is 1.69 × 10 m s , TOUGH2-MINC result can well fit the field BTC. Furthermore, the importance of matrix diffusion was verified by fitting the field BTC with analytical solutions that either neglect matrix diffusion or consider the mass exchange between mobile and immobile zones within the fracture as the attenuation transport mechanism. It was found that the BTC tailing can only be modeled if matrix diffusion is taken into account. Therefore, matrix diffusion can explain BTC tailing for advection-dominated tracer transport in fractured sandstone. (

(D) is 1.69 × 10 -6 m 2 s -1 , TOUGH2-MINC result can well fit the field BTC.Furthermore, the importance of matrix diffusion was verified by fitting the field BTC with analytical solutions that either neglect matrix diffusion or consider the mass exchange between mobile and immobile zones within the fracture as the attenuation transport mechanism.It was found that the BTC tailing can only be modeled if matrix diffusion is taken into account.Therefore, matrix diffusion can explain BTC tailing for advection-dominated tracer transport in fractured sandstone.

INTRODUCTION
Fluid flow and solute transport in fractured rock have attracted a lot of attention in recent years.These subjects are closely related to various fields of application, e.g., energy extraction in geothermal reservoirs, the fate and transport of contaminants in the subsurface, and the final disposal of spent nuclear fuels in fractured media.The last application concerns, among others, how radionuclides migrate through the interconnected fracture network.This problem can be analyzed by numerical modeling of tracer transport through fractured media.However, the modeling technique is challenged by the characterization of heterogeneity within the fracture network.In other words, a successful modeling of tracer transport in fractured media relies on a realizable description of the inherent heterogeneity and associated uncertainties of the complex, interconnected fracture network.
Despite the complexity of naturally fractured media, field evidence has shown that natural fracture networks can often be statistically classified into a finite number of fracture sets, e.g., Sanford et al. (2002).Accordingly, the entire network can be approximated as a few sets of fracture systems with rock matrix imbedded in between.For short-term (at time scales no longer than a few weeks) tracer experiments in such systems, tracer may still enter the rock matrix by molecular diffusion but is not expected to penetrate deeply enough to be influenced by adjacent fractures.Under this condition, the rock matrix can be treated as an infinite medium relative to a single fracture.In other words, the parallel fracture system can be further simplified as a semi-infinite system with a single fracture and an infinite rock matrix.This simplified conceptualization has been successfully applied for interpreting in-situ tracer tests in fractured media (Maloszewski andZuber 1990, 1993;Himmelsbach et al. 1998).
A typical phenomenon of tracer experiments in fractured media is the pronounced tailing of the breakthrough curve (BTC).This tailing implies the existence of a certain attenuation mechanism.Various mechanisms have been proposed to explain BTC tailing.Matrix diffusion has long been hypothesized as the main attenuation mechanism (Neretnieks 1980;Maloszewski and Zuber 1990;Moench 1995;Novakowski et al. 1995;Sanford et al. 1996).Other diffusive mechanisms similar to matrix diffusion have also been proposed, e.g., diffusion into interchannel small-aperture spaces in fractures (Abelin et al. 1994); advection-dispersion model with tran-sient solute storage (hereafter referred to as ADTS) conceptualizing mass exchange between mobile fluid in fractures and stagnant fluid along irregular fracture walls (Raven et al. 1988); and diffusion into an unspecified distribution of stagnant fluids associated with the fracture and matrix (Haggerty et al. 1999).On the contrary, Becker and Shapiro (2000) have verified that diffusive mechanisms fail to interpret their dipole and convergent tracer tests conducted in fractured crystalline rock that is characterized by a highly heterogeneous hydraulic conductivity field.They found that advective transport processes should be the only nondiffusive mechanism for the observed BTC tailing.These evidences indicate that modeling of solute transport in fractured rock is further complicated by conceptualizing the attenuation mechanism that leads to BTC tailing.
Fractured media are featured by their distinct hydraulic as well as transport properties between the fracture and rock matrix.That is, the fracture has high hydraulic conductivity but low porosity, whereas the matrix has low hydraulic conductivity but large porosity.Without loss of precision, it can be assumed that groundwater flow occurs entirely within the interconnected fracture network.Consequently, solute travels quickly through fractures by advection; whereas it migrates at a much slower rate into and within the rock matrix by molecular diffusion.Concentration in the matrix therefore distributes in a direction normal to the fracture-matrix interface as a decreasing function of the distance from the fracture.Taking the average of the concentration in the ensemble matrix apparently underestimates the concentration near the fracture-matrix interface.Unfortunately, this ensemble average is the concept of the quasi-steady state assumption adopted in the conventional double porosity model (DPM).Therefore, the DPM is limited in its accuracy of predicting the interporosity flux.For advection dominated transport cases, this underestimation tends to reduce the BTC tailing and leads to an unrealistic conclusion that matrix diffusion is insignificant.In order to accurately estimate interporosity flux, this paper uses the multiple interacting continua (MINC) (Wu and Pruess 2000) concept to enhance the resolution of spatial discretization near the fracture-matrix interface.By doing so, we not only get a reasonable estimation of concentration in the matrix but also increase the accuracy of predicting interporosity flux.
Besides the accuracy of simulating matrix diffusion, another difficulty of finite difference modeling of advection dominated transport is the numerical dispersion problem.Peaceman (1977) showed that numerical dispersion of finite difference modeling of the advection equation is a function of grid block size and time step size.Therefore, the advection-dispersion equation is mathematically equivalent to the advection equation if numerical dispersion can be calibrated to the same value as the hydrodynamic dispersion coefficient in the advection-dispersion equation.Migration of tracers in fractures is driven by advection as well as hydrodynamic dispersion, which can be described by the advection-dispersion equation.Therefore, successful finite-difference modeling of advection dominated tracer transport in fractured media relies on appropriate spatial as well as temporal discretizations.
This paper considers a short-term and radially convergent tracer test conducted in fractured sandstone.Unlike most tracer experiments that are commonly featured by moderate advection effects and marked BTC tailing, this dataset exhibits a strong advection effect yet a non-negligible BTC tailing.Since in-situ interference test results indicate that the fracture network seems to be dominated by a single fracture, it is expected that the hydraulic conductivity field is not very heterogeneous.This implies that a nondiffusive mechanism may not be important in the observed tailing.Hence, the objectives of this paper are to demonstrate the following perspectives through numerical simulation: (a) BTC tailing of in-situ tracer experimentation can be explained by matrix diffusion, even when the present experiment is apparently dominated by advection effect.(b) Accurate simulation of matrix diffusion can be done numerically by adopting a transient interporosity flux concept.This concept is implemented by incorporating the MINC discretization algorithm with the transport simulator TOUGH2 (Pruess et al. 1999).MINC is a refinement of the conventional DPM in that the dependence of tracer concentration in the matrix on a distance that is normal to the fracture-matrix interface is taken into account.
That is, a more realistic estimation of interporosity flux can be obtained, instead of the "averaged" flux estimated by DPM.(c) Although other attenuation mechanisms may have been proposed in the literature to explain BTC tailing, this paper will demonstrate that matrix diffusion is more likely to be the one that controls BTC tailing observed in the field.(d) The problem of numerical dispersion can be successfully overcome by suitably adjusting the time step size and grid block size such that simulation of the effect of hydrodynamic dispersion can be replaced by controlling the magnitude of numerical dispersion.

TRACER EXPERIMENT
The study site of this tracer experiment is located at a diversion tunnel of a power plant in Hoping village of Taichung County, Taiwan, which is about 5 kilometers west to Kukuan. Figure 1 shows the location of this test site, in which Taichung County is illustrated as the light grey area.
Rock formation at this site belongs to the Lileng member of the Paileng formation that is distributed from the north bank of the Tachiahsi basin to the west of Kukuan, see the dark grey area in Fig. 1.The Paileng formation consists of thick to massive, lower grade metamorphosed quartzitic sandstone that is grayish white or white in color.The thickness of the Lileng member is about 400 to 600 m.In a petrographic sense, the Paileng formation is equivalent to the Szeleng sandstone formation (CGS 1999).
Two parallel and inclined open boreholes (BH-2 and BH-3) with trend and plunge at 82°a nd 40°, respectively, were drilled on the eastern wall of the diversion tunnel.They are separated 5 m horizontally, and are both 30 m long in the dipping direction.See the upper inset in Fig. 1 for the relative position of the two boreholes.Totally seven packed-off sections were selected in these boreholes for subsequent hydraulic and tracer tests.Depth intervals of these sections are listed in Table 1.During interference and tracer tests, a multiple packer system was installed in BH-2 for eliminating hydraulic interaction between packed-off sections, and a double packer system was installed in BH-3 at a particular packed-off section.Pumping test results showed that hydraulic conductivities at the upper two packed-off sections of BH-3 are approximately 10 -5 and 10 -4 m s -1 for the bottom section of BH-3.
Both long-term monitored pressure data and cross-hole interference test results indicate a noticeable interconnectedness of the fracture network between BH-2 and BH-3, especially at a depth of approximately 25 m.This is supported by the following observations: First of all, despite temporal fluctuations of monitored pressure data, all packed-off sections in BH-2 show response shortly after rainfall events.This is revealed by an elevated pressure increase in all pressure-time curves.Comparisons with precipitation data collected at a gauge station near the test site have verified the link between pressure increase and precipitation.Hence, interconnectedness of the fracture network is likely to extend from the ground surface to as deep as 25 m.Secondly, if BH-3 is open and pumped continuously at a constant rate, all packed-off sections in BH-2 respond accordingly to the pumping, though at an irregular frequency.In particular, the packed-off section at 25 m in BH-2 shows the largest pressure increase.Thirdly, if BH-3 is Table 1.Depths of packed-off sections in boreholes.
packed-off and pumped at any one of the three sections, the lower two sections in BH-2 respond to the pumping more significantly than the upper two sections.Again, the bottom section in BH-2 shows the largest pressure increase.Accordingly, the rock body between BH-2 and BH-3 is better connected in the lower part than in the upper part.
Although interference test results cannot provide a detailed network structure between BH-2 and BH-3, these data do confirm the fact that BH-2 is well connected to BH-3 through some water conducting features (WCFs) at a depth of 25 m.Here, WCF is defined as a hydrogeological structure with enhanced transmissivity within a rock body (Mazurek 2000).Moreover, interference test results suggest that the rock mass at this depth may be dominated by a parallel fracture system (Chiang et al. 2001).Accordingly, it was decided to inject tracer into the lowermost section of BH-2 and withdraw water from the lowermost section of BH-3, both of which are approximately at a depth of 25 m.
A conservative tracer, sodium chloride (NaCl), was used in this experiment.It was prepared by diluting 1 kg of NaCl in 20 liters of water, and was injected as a pulse into BH-2.The pumping rate was fixed at 1.5 l min -1 .Before starting the tracer test, a steady state flow field was created in the field by pumping BH-3 continuously until pressure in BH-3 became stabilized.Electrical conductivity of the withdrawn solution was manually measured in the field, from which tracer concentration could be back-calculated by the following linear calibrating equation (Chiang et al. 2001) in which C (kg m -3 ) and EC ( µmho cm -1 ) are the concentration and electrical conductivity of the NaCl solution, respectively.Other mathematical symbols used in this paper are defined in the 'Notation' section in the Appendix.Tracer BTC obtained from the field experiment is shown in Fig. 2. The initial plateau of the BTC reflects the electrical conductivity of local groundwater, which is approximately 180 µmho cm -1 .This value is considered to be the background electrical conductivity of the in-situ groundwater.It is added to the injected tracer concentration to yield the tracer BTC in Fig. 2. The first arrival, the peak, and the start of tailing occur at approximately 300, 340, and 390 minutes after injection, respectively.Qualitatively, this means that tracer transport at this site is dominated by advection.Note that attenuated tracer transport is indicated by the slight asymmetry in the late-time BTC.In summary, the present tracer experiment is characterized by a strong advection effect in the fracture and delayed transport that is caused by a specific attenuation mechanism.

Conceptual Model
Interference test results imply that a quick flow path formed by parallel fractures may exist between BH-2 and BH-3.Whether the interaction between adjacent flow paths is significant enough to influence tracer transport in an individual flow path is a concern.Another issue is that whether a unidimensional coordinate system is adequate to approximate the radial convergent tracer experiment.First of all, note that this experiment lasted only for a few hours.Within such a short period of time, tracer may not have enough time to migrate into the tight rock matrix deeply enough to be influenced by adjacent fractures.Secondly, tracer transport at this site was dominated by advection.That is, transverse dispersion is expected to be negligible compared with advective transport.Although interference test results and pressure monitoring data show vertical connection within the rock body, its effect on tracer transport is not discussed in this paper.Therefore, a unidimensional coordinate system with a single fracture and a semi-infinite rock matrix is conceptualized as an approximation of the test conditions.Maloszewski andZuber (1990, 1993) have validated this single fracture conceptual model by successfully interpreting several short-term to intermediate-term tracer test data with their single fissure dispersion model (SFDM) solution.Moreover, Himmelsbach et al. (1998) have confirmed the applicability of SFDM to short-term tracer experiments in fractured granite at the Lindau test site located in the southern Black Forest of southwestern Germany.Note that simplifying the radial flow system as a unidimensional one means that the radial flow field can be lumped into a unidimensional flow field in a dominant single fracture.Also, parameters obtained from a single fracture model thus represent lumped averages of fracture parameters of the rock body under investigation.
The conceptual model of this single fracture system is illustrated in Fig. 3.A unidimensional single fracture with smooth surfaces separated by a constant aperture of 2b is considered, which is sandwiched between semi-infinite rock matrix.A steady state groundwater flow with a constant pore velocity (V) is assumed in the fracture, whereas no groundwater flow occurs in the matrix media.Because surface roughness of the fracture is neglected and fracture aperture is generally small (on a millimeter scale or smaller) for sandstones, complete mixing of tracer across the entire width of the fracture can be assumed.Tracer transport along the fracture is assumed to be governed by advection as well as hydrodynamic dispersion.Furthermore, tracer transport across the fracture-matrix interface and within the granular pore space of the rock matrix is assumed to be controlled by molecular diffusion only.In Fig. 3, rock matrix affected by matrix diffusion is schematically illustrated by the dotted region, and curved arrows represent molecular diffusion in the porous matrix.

Numerical Scheme
Numerical simulations of tracer transport in a single fracture system were carried out by the general-purpose TOUGH2 code, a numerical simulator for non-isothermal, multiphase and multi-component fluid flow and solute transport in one, two, and three-dimensional porous and fractured media (Pruess et al. 1999).
Unlike most simulators that solve mass balance equations in a differential form, TOUGH2 solves the following integral mass balance equation: in which the superscript κ is the component index; M ( ) κ is the mass of component κ per unit volume of grid block V n ; F ( )  κ is the sum of the κ -th component advective and dispersive fluxes that are normal to the interface Γ n , and q ( ) κ is the sink/source term (including adsorption as well as diffusive mass loss at the fracture-matrix interface) of component κ per unit volume of grid block.Equation ( 2) can be expanded into NK equations if NK fluid components are considered.For isothermal flow simulations, TOUOGH2 solves Eq. ( 2) by recasting M ( ) κ and F ( ) κ in terms of fluid saturation and pore pressure, respectively.For isothermal simulations, Fig. 3.The single fracture conceptual model.

M ( )
κ depends additionally on the mass fraction of the κ -th fluid component.Fluid saturation, pore pressure, and component mass fraction are therefore called primary variables in TOUGH2 which are solved by a user-selected linear equation solution package imbedded in TOUGH2.Since isothermal, single aqueous phase (water) and two-component (water and tracer) solute transport in fully saturated fractured media is considered in this paper, the primary variables are pore pressure and component mass fraction.Initial and boundary conditions can then be defined in terms of these two primary variables.
Using divergence theorem, the surface integral in Eq. ( 2) can be transformed into a volume integral.Considering a one-dimensional uniform flow field, the TOUGH2 governing equation can be expressed in a form similar to the conventional advection-diffusion equation: where C ( ) κ , V, and D are the concentration of component κ , steady pore velocity, and hydro- dynamic dispersion coefficient, respectively.In deriving Eq. ( 3), it is assumed that Darcy's law is valid and that dispersive as well as diffusive fluxes follow Fick's first law.Adsorption is assumed to follow a linear equilibrium isotherm.If only adsorption is considered in the sink/ source term, then Eq. ( 3) can be further simplified by dropping the term q ( ) κ and replacing , in which R is the retardation factor.For conservative tracers such as NaCl, the retardation factor is reduced to one.
Integral finite difference method is employed in TOUGH2 for spatial discretization, while fully implicit time stepping is adopted for achieving unconditional stability (Pruess et al. 1999).The dispersive mass flux in TOUGH2 can be simulated by the two-dimensional transport module known as T2DM (Oldenburg 1993), in which the unknown Darcy velocities parallel to grid block boundaries are interpolated from known Darcy velocities by the 'distance-averaged' scheme in order to calculate the hydrodynamic dispersion tensor.However, simulation of hydrodynamic dispersion for the present advection-dominated case is replaced by a numerical skill, which will be explained in section 4.2.
The conceptual model in Fig. 3 assumes that matrix diffusion is the sole attenuation mechanism.When advection dominates tracer transport in the system, it is anticipated that the accuracy of calculating concentration flux across the fracture-matrix interface (so-called interporosity flux) must be high enough to simulate the matrix diffusion term.The conventional DPM may not be feasible because it tends to underestimate this interporosity flux due to its quasi-steady-state assumption.In this paper, therefore, the MINC (Pruess et al. 1999) concept is adopted.The fundamental concept of MINC is based on the observation that changes of concentration in fractured media are expected to propagate more rapidly through the fracture, while only slowly invading the matrix (Wu and Pruess 2000).Accordingly, distribution of tracer concentration in the matrix will be controlled locally by the perpendicular distance between the matrix and fracture.In a numerical sense, MINC is equivalent to discretizing the rock matrix into several continua according to the perpendicular distance from the nearest fracture element.In MINC, the discretization is done numerically by providing volumetric fractions of the matrix.Arbitrary volume fractions of the matrix continua can be given.However, it is better to specify gradually increasing fractions away from the fracture such that the spatial variability of interporosity flux can be better resolved.Since the subgridding step in MINC follows the same format as the integral finite difference method, MINC can be directly incorporated into TOUGH2 without extra code changes.A schematic example of MINC subgridding is illustrated in Fig. 3 by the dashed lines, in which four matrix continua are additionally discretized at either side of the fracture continuum.Note that MINC reduces to DPM if the rock matrix is lumped into a single continuum.Consequently, the interporosity flux in DPM is proportional to the difference of lumped averages of tracer concentrations in the fractured media and in rock matrix.

RESULTS AND DISCUSSIONS
In order to clarify the approach adopted in this paper, the key simulation procedure is summarized in the flowchart shown in Fig. 4. Detailed explanations of each individual element can be seen in the following sections.

Initial and Boundary
Numerical grid of all TOUGH2-MINC simulations is designed by first generating a onedimensional grid with a constant grid size as if there is only one single porous continuum.This continuum is then further discretized in a direction transverse to the flow direction by MINC according to specified volume fractions.The first fraction corresponds to fracture porosity, and the rest of the fractions designate how the matrix continuum is partitioned.In order to enhance the resolution of interporosity flux at the fracture-matrix interface, volume fractions of the matrix are gradually increased away from the fracture.Based on parsimony principles, totally six matrix continua are discretized because simulation results become stabilized if more than six matrix continua are considered.
All simulations start from obtaining the uniform flow field and then simulating tracer transport in the system.First, upstream and downstream pressures are calculated according to the given pumping rate.The uniform flow field is then created by running TOUGH2-MINC to a steady-state.Next, the steady-state flow field is substituted into transport simulations, with a zero initial concentration throughout the system.The pulse injection is simulated by introducing an instantaneous tracer mass at the grid block next to the upstream grid block.Hydrodynamic dispersion can be simulated by incorporating the T2DM module (Oldenburg 1993) into TOUGH2-MINC.Because of the numerical dispersion problem, however, this mechanism is simulated by the step size refinement method discussed in the next section.Finally, background electrical conductivity is converted to concentration according to Eq. ( 1) and is added to simulated tracer concentration to yield the tracer BTC.

Numerical Dispersion
Note that finite-difference simulation of advection-dominated transport inevitably deteriorates due to numerical dispersion.Usually, the significance of advection is quantified by the magnitude of the Peclet number (Pe) which is defined as: in which L is the interwell distance and α L is the longitudinal dispersivity.The larger the advection effect, the greater the Pe.Peaceman (1977) has derived the magnitude of numerical dispersion in finite-difference simulations for various spatial and temporal discretization schemes.Since TOUGH2 uses an upstream weighting method in space and fully implicit discretization in time, numerical dispersion (D num ) should be V x ∆ ( )/ 1 2 + λ , where λ is the Courant number that is defined as V x ∆ ∆ t / in which ∆t and ∆x are time step size and grid size, respectively.Because V and ∆x are fixed, D num thus increases with ∆t .
In order to simulate hydrodynamic dispersion with the presence of numerical dispersion, TOUGH2-MINC simulations were carried out by turning off the physical hydrodynamic dispersion (letting α L = 0 ) and refining the time step size until the magnitude of numerical disper- sion is equal to the expected hydrodynamic dispersion coefficient.In other words, numerical simulation of the advection-dispersion equation, Eq. ( 3), is reduced to simulation of the math-ematically equivalent advection equation by controlling the magnitude of numerical dispersion.
To validate the refinement method, first note that the MINC concept reduces to a single porosity model, i.e., a porous medium.Then, TOUGH2-MINC simulations can be treated as simulating tracer transport in porous media.Lenda and Zuber (1970) proposed a dispersion model (hereafter referred to as DM) that considers tracer transport in porous media with a pulse injection.Therefore, this refinement method can be validated by comparing the TOUGH2-MINC result with the DM solution that uses the same hydrodynamic dispersion coefficient as the numerical dispersion in TOUGH2.Results show that this refinement method works well despite the magnitude of Pe.

Simulation Results
Average values of fracture aperture and matrix porosity of the Szeleng sandstone formation were investigated by Cheng (2001) and Lin (2002) as 200 µm and 2%, respectively.These values were used as input data for the base case of TOUGH2-MINC simulations.Sensitivity of simulation results on these parameters will be studied later in this section.
The steady flow field in the single-fracture system is induced by pumping.The resultant constant pore velocity is estimated by dividing the interwell distance by the peak arrival time of tracer (t 0 ).From Darcy's law, the relationship between pore velocity and pore pressure difference ( ∆P ) is shown in Eq. ( 5), in which k f , µ , and n f are fracture permeability, water viscosity and fracture porosity, respectively.
Note that the gravity term ( ρg ) in Eq. ( 5) is neglected because a single fracture is assumed in the horizontal direction.Thus, pore pressure difference can be calculated by Eq. ( 5), which is then supplied to TOUGH2 as the pressure boundary condition.Input data for the base case of TOUGH2-MINC simulations are summarized in Table 2. Injected tracer mass is estimated from the recovery data.From the field BTC, tracer recovery can be estimated as: Substituting field BTC data into Eq.( 6) yields the recovered tracer mass to be approximately 14.4 g.Since interference test results and pressure monitoring data show vertical connection within the rock body, tracer mass may deviate from the single fracture system and results in incomplete recovery in the field.Downstream pore pressure (P d ) considers hydrostatic pressure at a depth of 25 m, which is added to the pressure difference calculated from Eq. ( 5) to give the upstream pore pressure (P u ).Fracture permeability is converted from the hydraulic conductivity estimated from the field pumping test result, which is approximately 10 -11 m 2 .Its relationship with fracture porosity is by Zuber's model (Zuber 1974), in which a network with tortuous fissures of the same width is considered, i.e.,: where τ f is the tortuosity factor in the fracture.Since k f and fracture aperture are 10 -11 m 2 and 200 µm, respectively, fracture porosity can be calculated from Eq. ( 7) as 0.3% if τ f is assumed to be one.Values of time step size ( ∆t ) and molecular diffusion in matrix (D m ) are selected such that the simulated BTC is close to the field BTC.Using the input data in Table 2, Fig. 5 compares simulated BTC and the field BTC.Obviously, the simulated BTC fits the field data well.
Using the same parameters as in Table 2 but varying D m and ∆t , Fig. 6 shows the sensitivity of each parameter on simulation results.As D m decreases, the effect of matrix diffusion diminishes and more tracer mass remains in the fracture.Hence, compared with the base case BTC, the resultant BTC (the dash-dot-dot curve in Fig. 6a) has a higher peak concentration.In addition, this BTC is almost symmetric.That is, the symmetry of a BTC reflects the declining effect of matrix diffusion in the fractured rock.Conversely, increasing D m will enhance matrix diffusion Table 2. Input data of the base case of TOUGH2-MINC simulations.and result in pronounced tailing, see the dashed curve in Fig. 6a.Therefore, matrix diffusion is not a negligible mechanism for tracer transport at this site and should be correlated with the observed BTC tailing.By the step size refinement method, simulation of hydrodynamic dispersion is replaced by adjusting ∆x and/or ∆t until the magnitude of numerical dispersion is equal to the expected hydrodynamic dispersion coefficient.It is already known that the magnitude of numerical dispersion is .Since ∆x is fixed in all TOUGH2-MINC simulations, D num thus increases with ∆t .Figure 6b compares three simu- lation results by increasing ∆t from 3 to 30 and 300 sec, corresponding to D = 8.44 × 10 -7 , 1.69 × 10 -6 , and 1.01 × 10 -5 m 2 s -1 , respectively.Changing hydrodynamic dispersion affects only the spreading of the BTC but not BTC tailing.Therefore, the three BTCs in Fig. 6b have approximately the same peak arrival time but the spreading increases as ∆t increases.Note that the base case ∆t (30 sec) corresponds to D as 1.69 × 10 -6 m 2 s -1 .Substituting this value into Eq.( 4) yields the Pe as 740, a very large value manifesting the dominance of advection in the field.
Fracture aperture (2b) and matrix porosity (n m ) have opposite effects on tracer transport.Figure 7 compares simulation results using different values of 2b and n m , with the same steady flow field in all cases.By increasing 2b or reducing n m , more tracer mass will be preserved in the fracture than is diffused into the rock matrix.Hence, the corresponding BTC has a higher peak concentration and less tailing than the base case BTC, see the dash-dot-dot curve and delta symbol in Fig. 7. On the contrary, reducing 2b or increasing n m allows more tracer mass to be diffused into the porous matrix.This leads to a more elevated tailing effect than the base case BTC, see the dashed curve and square symbol in Fig. 7. Maloszewski and Zuber (1993) combine n m , D m and 2b into a non-dimensionless "diffusion parameter" to consider the effect of matrix diffusion for conservative tracers.This diffusion parameter (a) is formulated as: Physical meaning of diffusion parameter is that the larger the diffusion parameter, the greater the matrix diffusion.Therefore, increasing n m or decreasing 2b will enhance matrix diffusion, and vice versa.From Eq. ( 8), the diffusion parameter of the base case (a 0 ) is 3.16 × 10 -4 sec -0.5 .Values of 2b and n m in Fig. 7 are selected such that the new diffusion parameter is either 4 or 0.25 a 0 .If matrix diffusion becomes stronger, i.e., a is increased to 4 a 0 , Fig. 7 shows that the corresponding BTCs have a longer tail.On the contrary, BTCs in Fig. 7 corresponding to 0.25 a 0 are more symmetric and have a larger peak concentration, meaning a less significant effect of matrix diffusion.Therefore, the mechanism that controls the tailing or peaking of the BTC can be ascribed to matrix diffusion.

Comparisons with Other Models
In the following, simulation results from TOUGH2-MINC are compared with those from DM, DPM, and ADTS.The comparisons aim at clarifying the following facts that (1) rock matrix plays the role of delaying tracer transport in the fracture because tracer mass may enter the matrix by a much slower mechanism of molecular diffusion, (2) DPM lacks the accuracy of estimating interporosity flux because of the quasi-steady state assumption, and (3) matrix diffusion is better than the mass exchange mechanism in ADTS to explain BTC tailing.
In order to quantitatively compare fitting results by different models, the root-mean-squareerror (RMSE) between field BTC and simulated BTC is computed, which is defined as: Note that C Field in Eq. ( 9) is interpolated from the simulated BTC such that it corresponds to the same time at which the field BTC is observed.In addition, the RMSE corresponding to times greater than 5.5 hours is also computed in order to emphasize the model fit at the BTC tail.
Computed RMSEs for different models are listed in Table 3.

DM
It is clear from Fig. 6 that BTC tailing is reasonably captured by TOUGH2-MINC because it considers attenuated transport by matrix diffusion.To demonstrate the effect of matrix diffusion, the in-situ BTC was also fitted with the DM solution.Figure 8a (see the curve labeled with the ' ×' symbol) shows that the resultant BTC cannot fit the field BTC because it is perfectly symmetrical and shows no tailing at all.This is because DM does not take matrix diffusion into account.The poor fit of DM is also demonstrated by its RMSE.Note that DM has the largest RMSE compared with all other models.
Theoretically, TOUGH2-MINC reduces to DM if the diffusion parameter is zero.By letting n m = 0 or D m = 0, a perfectly symmetric BTC is obtained from TOUGH2-MINC.This  Table 3. Root mean squared errors for all simulation results.
BTC is practically the same as the one modeled by DM.This result again verifies that BTC tailing cannot be simulated if the attenuation mechanism in terms of matrix diffusion is neglected.
In other words, the importance of matrix diffusion is justified by numerical simulation.

DPM
Note that MINC reduces to DPM if only one matrix continuum is considered.Since DPM uses the quasi-steady state assumption and considers the rock matrix as a single continuum that is relatively infinite to the fracture continuum, it tends to underestimate interporosity flux.For example, the DPM BTC in Fig. 8a (see the dashed curve) practically shows no tailing and is very close to the DM BTC.Since DM completely neglects matrix diffusion, this result reveals that DPM cannot properly simulate the effect of matrix diffusion if advective transport is important.Also, RMSEs calculated from DPM are all greater than those from TOUGH2-MINC, indicating again the inadequacy of DPM in simulating the interporosity flux.
If interporosity flux calculated by MINC (f MINC ) and DPM (f DPM ) are compared, it was found that the ratio f f MINC DPM can be as large as 10 9 .Underestimation of interporosity flux further affects simulation of tracer transport in the rock matrix.Although no field measurements of tracer concentration in the matrix were available, simulation results in Fig. 8b predict that DPM BTC in the rock matrix is at least four orders of magnitude smaller than the one simulated by TOUH2-MINC.For simplicity, only concentration values greater than 10 -8 kg m -3 are shown in Fig. 8b.Also, Fig. 8b shows that matrix concentration increases quickly to the peak value and then declines at a much smaller rate after peak concentration.Note that the peak concentration in the matrix occurs shortly after the peak concentration in fracture.In addition, the declination persists even after the tracer is completely removed from the fracture.This result therefore supports the idea by Sudicky and Frind (1982) that the matrix block plays the role of dynamic mass storage that constitutes an attenuation mechanism for the advance of tracer in the fracture system.

ADTS
It is mentioned in the introduction that Raven et al. (1988) used the advection-dispersion model with transient solute storage in immobile fluid zones to interpret the pronounced BTC skewness observed in their field data.In order to investigate whether transient solute storage is an effective attenuation mechanism for explaining the present field data, the in-situ BTC was also fitted with the ADTS solution.The optimal ADTS BTC shown in Fig. 8a (see the dashdot-dot curve) uses the following fitting parameters: Pe = 360, φ = 0.95, N = 60, and τ = 5.55 hrs, in whichφ , N, and τ are the fraction of mobile fluid phase to the total fracture pore space, dimensionless mass exchange constant between mobile and immobile fluid phases, and plug flow transient time, respectively.The RMSEs corresponding to this BTC are shown in the last row of Table 3.Although ADTS fits the tail better than TOUGH2-MINC, its overall RMSE is still greater than that by TOUGH2-MINC.That is, the rising limb of the field BTC cannot be fitted by ADTS.
The essence of ADTS is the existence of an immobile fluid phase, from which solute mass can be exchanged between mobile and immobile fluid phases and thus results in BTC skewness.Hence, the smaller the φ value, the larger is the BTC skewness.However, the optimal ADTS BTC uses a large φ value.This means that the fracture pore space is almost filled with the mobile fluid phase.In fact, if the maximum φ was used, i.e., φ = 1.0, the resultant BTC is practically the same as the one shown in Fig. 8a.On the contrary, if φ is smaller than 0.9, the resultant BTC becomes too dispersed to fit the field BTC.That the optimal ADTS uses a large φ value is not a surprise because ADTS is effective in accounting for pronounced BTC skew- ness due to attenuated transport by the presence of the immobile fluid phase.Since the field data does not have a strong BTC skewness, it is anticipated that ADTS should use a large φ value.In other words, the transient solute storage concept is not appropriate to explain the field BTC.

CONCLUSIONS
A short-term tracer experiment conducted in fractured sandstone was analyzed.Tracer transport at the site was obviously dominated by advection.Yet, a non-negligible BTC tailing was also observed.Based on field interference test results, a single fracture conceptual model was proposed.In this model, rock matrix is considered to be infinite relative to the single fracture because of the short duration of the testing.In addition, tracer is assumed to travel through the single fracture by advection as well as hydrodynamic dispersion.Molecular diffusion is assumed to be the attenuation mechanism that controls the migration of tracer between the fracture and rock matrix.Because of overwhelming advection in the fracture, interporosity flux across the fracture-matrix interface has to be accurately calculated in order to reflect the observed BTC tailing.The MINC concept has proved to be effective in modeling this matrix diffusion.In fact, modeling results predict that the interporosity flux simulated by DPM is order-of-magnitudely smaller than that by TOUGH2-MINC.
The field BTC was compared with the DM and DPM solutions.DM does not consider matrix diffusion.Although DPM considers matrix diffusion, its quasi-steady state assumption inevitably underestimates the interporosity flux.It was found that the DM result is practically the same as the DPM result, both models are less efficient than TOUGH2-MINC to simulate the observed BTC tailing.Moreover, an alternative attenuation mechanism, ADTS, considering transient solute storage was also considered.Results showed that the ADTS BTC is almost symmetric because of the large fraction of mobile fluid to the total fracture pore space.Consequently, transient solute storage due to immobile fluid phase is not appropriate to explain the field data.Hence, the assumption that matrix diffusion is the key attenuation mechanism for the observed BTC tailing is verified.
Despite the simplicity of the conceptual model, matrix diffusion as the attenuation mechanism can be proved to be effective in explaining the field data.Although estimated fracture aperture and matrix porosity measured from the Szeleng formation were used in the simulation, the good match between simulated results and field data suggests that averaged transport parameters of the Szeleng formation should be representative of the transport characteristics of the Paileng formation.Liou, T. S., 2007: Numerical analysis of a short-term tracer experiment in fractured sandstone.

Fig. 1 .
Fig. 1.Location map of the study site.

Fig. 6 .
Fig. 6.Sensitivity of TOUGH2-MINC with respect to (a) molecular diffusion coefficient in matrix and (b) time step size.

Fig. 8 .
Fig. 8.Comparison of TOUGH2-MINC and other transport models: (a) simulated and field BTCs in the fracture and (b) simulated BTCs in the rock matrix at the withdrawal borehole.

C
Field tracer concentration of the field BTC (M L -3 ) C' interpolated tracer concentration from simulated BTC (M L -3 ) D hydrodynamic dispersion coefficient in fracture (L 2 T -1 ) D m molecular diffusion coefficient in rock matrix (L 2 T -1 ) D num numerical dispersion (L 2 T -1 ) EC electrical conductivity of tracer ( µmho L -1 ) F ( ) κ sum of interfacial mass fluxes of component κ (M L -2 T -1 ) component fluid mass per unit volume of grid block (M L -3 ) N dimensionless mass exchange constant between mobile and immobile fluid phases P u upstream pore pressure (M L -1 T -2 ) P d downstream pore pressure (M L -1 T -2 ) Pe Peclet number Q steady volumetric flow rate through the fracture (L 3 T -1 ) R retardation factor V uniform pore velocity in fracture (L T -1 ) V volume of a general grid block (L 3 ) V n volume of the n-th grid block (L 3 ) a diffusion parameter (T -0term of component κ per unit volume of grid block (M L -3 ) t time (T) t 0 mean transit time of water (T) α L longitudinal dispersivity (L) ∆x grid block size (L) ∆t time step size (T) ∆P pressure difference between upstream and downstream ends (M L -1 T -2 ) φ fraction of mobile fluid phase to the total fracture pore space λ general grid block (L 2 ) Γ n interface of the n-th grid block (L 2 )