An interpretation of potential scale dependence of the effective matrix diffusion coefficient

Matrix diffusion is an important process for solute transport in fractured rock, and the matrix diffusion coefficient is a key parameter for describing this process. Previous studies have indicated that the effective matrix diffusion coefficient values, obtained from a large number of field tracer tests, are enhanced in comparison with local values and may increase with test scale. In this study, we have performed numerical experiments to investigate potential mechanisms behind possible scale-dependent behavior. The focus of the experiments is on solute transport in flow paths having geometries consistent with percolation theories and characterized by multiple local flow loops formed mainly by small-scale fractures. The water velocity distribution through a flow path was determined using discrete fracture network flow simulations, and solute transport was calculated using a previously derived impulse-response function and a particle-tracking scheme. Values for effective (or up-scaled) transport parameters were obtained by matching breakthrough curves from numerical experiments with an analytical solution for solute transport along a single fracture. Results indicate that a combination of local flow loops and the associated matrix diffusion process, together with scaling properties in flow path geometry, seems to be an important mechanism causing the observed scale dependence of the effective matrix diffusion coefficient (at a range of scales).     

Interpretation of injection-withdrawal tracer experiments conducted between two wells in a large single fracture

Tracer experiments conducted using a flow field established by injecting water into one borehole and withdrawing water from another are often used to establish connections and investigate dispersion in fractured rock. As a result of uncertainty in the uniqueness of existing models used for interpretation, this method has not been widely used to investigate more general transport processes including matrix diffusion or advective solute exchange between mobile and immobile zones of fluid. To explore the utility of the injection-withdrawal method as a general investigative tool and with the intent to resolve the transport processes in a discrete fracture, two tracer experiments were conducted using the injection-withdrawal configuration. The experiments were conducted in a fracture which has a large aperture (>500 microm) and horizontally pervades a dolostone formation. One experiment was conducted in the direction of the hydraulic gradient and the other in the direction opposite to the natural gradient. Two tracers having significantly different values of the free-water diffusion coefficient were used. To interpret the experiments, a hybrid numerical-analytical model was developed which accounts for the arcuate shape of the flow field, advection-dispersion in the fracture, diffusion into the matrix adjacent to the fracture, and the presence of natural flow in the fracture. The model was verified by comparison to a fully analytical solution and to a well-known finite-element model. Interpretation of the tracer experiments showed that when only one tracer, advection-dispersion, and matrix diffusion are considered, non-unique results were obtained. However, by using multiple tracers and by accounting for the presence of natural flow in the fracture, unique interpretations were obtained in which a single value of matrix porosity was estimated from the results of both experiments. The estimate of porosity agrees well with independent measurements of porosity obtained from core samples. This suggests that: (i) the injection-withdrawal method is a viable tool for the investigation of general transport processes provided all relevant experimental conditions are considered and multiple conservative tracers are used; and (ii) for the conditions of the experiments conducted in this study, the dominant mechanism for exchange of solute between the fracture and surrounding medium is matrix diffusion.     

Prediction of some in situ tracer tests with sorbing tracers using independent data

Some recent converging tracer tests with sorbing tracers at the Asp? Hard Rock Laboratory in Sweden, the TRUE tests, have been predicted using only laboratory data and hydraulic data from borehole measurements. No model parameters were adjusted to obtain a better fit with the experiments. The independent data were fracture frequency and transmissivity data obtained in the field and laboratory data on sorption and matrix diffusion. Transmissivity measurements in five boreholes in the rock volume containing the region surrounding the injection and collection points show that there is a high frequency of water conducting fractures. Of 162 packed off sections with 0.5 m packer distances, 112 were found to have a transmissivity above the detection limit. The specific flow-wetted surface (FWS) of the rock mass could be estimated from these data. The transmissivities were found to be reasonably well described by a lognormal distribution. Laboratory data on diffusion and sorption properties together with the hydraulic data were used to "predict" the residence time distribution (RTD) of the sorbing tracers. The results were compared with the experimental breakthrough curves. In these experiments, the water residence time is very small compared to the residence time of the sorbing tracers due to their diffusion and sorption within the rock matrix. We thus could neglect the influence of the water residence time in our predictions. Therefore, no information on water residence times or on "dispersion" was needed. The dispersion of the sorbing tracers is caused by the different sorbing tracer residence times in different pathways. The sorbing tracer residence time is determined by the ratio of flowrate to the flow-wetted surface in the different pathways and not by the water residence time. Assuming a three-dimensional flow pattern and using the observed fracture frequency and flowrate distribution, breakthrough curves for three strongly sorbing tracers were predicted. Only the laboratory data, the transmissivity measurements and the pumping flowrate were used in the predictions. No information on the water residence time as obtained by the nonsorbing tracers was used. The predictions were surprisingly accurate.     

Simulation and analysis of solute transport in 2D fracture/pipe networks: the SOLFRAC program

The Time Domain Random Walk (TDRW) method has been recently developed by Delay and Bodin [Delay, F. and Bodin, J., 2001. Time domain random walk method to simulate transport by advection-dispersion and matrix diffusion in fracture networks. Geophys. Res. Lett., 28(21): 4051-4054.] and Bodin et al. [Bodin, J., Porel, G. and Delay, F., 2003c. Simulation of solute transport in discrete fracture networks using the time domain random walk method. Earth Planet. Sci. Lett., 6566: 1-8.] for simulating solute transport in discrete fracture networks. It is assumed that the fracture network can reasonably be represented by a network of interconnected one-dimensional pipes (i.e. flow channels). Processes accounted for are: (1) advection and hydrodynamic dispersion in the channels, (2) matrix diffusion, (3) diffusion into stagnant zones within the fracture planes, (4) sorption reactions onto the fracture walls and in the matrix, (5) linear decay, and (6) mass sharing at fracture intersections. The TDRW method is handy and very efficient in terms of computation costs since it allows for the one-step calculation of the particle residence time in each bond of the network. This method has been programmed in C++, and efforts have been made to develop an efficient and user-friendly software, called SOLFRAC. This program is freely downloadable at the URL (labo.univ-poitiers.fr/hydrasa/intranet/telechargement.htm). It calculates solute transport into 2D pipe networks, while considering different types of injections and different concepts of local dispersion within each flow channel. Post-simulation analyses are also available, such as the mean velocity or the macroscopic dispersion at the scale of the entire network. The program may be used to evaluate how a given transport mechanism influences the macroscopic transport behaviour of fracture networks. It may also be used, as is the case, e.g., with analytical solutions, to interpret laboratory or field tracer test experiments performed in single fractures.     

Analysis of the migration of instantaneously injected cesium in artificial fractures of Lac du Bonnet granite, Manitoba, Canada

Breakthrough curves of ¹³⁷Cs and tritiated water injected instantaneously into artificial fractures in Lac du Bonnet granite were analyzed using the analytical solution for a single rock-fracture system and assuming the linear sorption isotherm of the solute. Parameters of nuclide diffusion and sorption in rock matrices, obtained by fitting, varied depending on the flow velocity in the fractures. According to theoretical calculations, different fracture flow velocities lead to different diffusion distances of nuclides in matrices at the same injection volume. As microscopic inhomogeneity is considered to exist in the rock matrix, the average diffusion-sorption characteristics of the matrix within the diffusion distance may have varied with the fracture flow velocity. Surface sorption was marked in fractures that had relatively high matrix sorption-diffusion capacities. The phenomenon was interpreted using the theoretical relationships developed between the surface sorption, matrix sorption and pore diffusion coefficient, and the porosity of matrices.The effect of the nonlinear sorption of solute was examined by numerically solving model equations that incorporate the nonlinear isotherm. This incorporation may contribute to the reduction of deviations between theoretical and experimental BTC's.     

Two-phase flow and transport in a single fracture-porous medium system

The prognosis for the remediation of contaminated fractured media is much worse than that for more homogeneous units. Fractures act as conduits for the flow of dense non-aqueous phase liquids (DNAPLs), while diffusion is responsible for the storage of dissolved mass in the surrounding matrix. A numerical model incorporating aqueous phase transport in a variable-aperture fracture and its surrounding matrix is developed and coupled with an existing two-phase flow model. The processes of transient two-phase flow, non-equilibrium dissolution, advective–dispersive transport in the fracture, and three-dimensional matrix diffusion are included in the model. Results from various investigations show that the DNAPL distribution is very sensitive to variations in aperture within a single fracture. Diffusion-controlled mass removal from both the matrix and from the hydraulically inaccessible zones within the fracture itself result in extremely large time frames for significant mass removal from these systems. Success in aqueous phase mass removal from the matrix is very sensitive to the effective fracture spacing. The hydraulic gradient in the fracture only affects the amount of water removed from the system, and does not greatly affect the amount of time required to remove the contaminant mass from the matrix. The ability to remove mass is somewhat sensitive to the porosity and effective matrix diffusion coefficient over the range of expected values.     

Radionuclide transport and retardation in rock fracture and crushed rock column experiments

Transport and retardation of non-sorbing tritiated water and chloride and slightly sorbing sodium was studied in Syyry area SY-KR7 mica gneiss, in altered porous tonalite and in fresh tonalite. Experiments were performed using dynamic fracture and crushed rock column methods. Static batch method for sodium was introduced to compare retardation values from static and dynamic experiments. The ¹⁴C-PMMA method was used to study the pore structure of matrices. The pore aperture distribution was evaluated from Hg-porosimetry determinations and the surface areas were determined using the B.E.T. method. The flow characteristics and transport behavior of tracers were interpreted using a numerical compartment model for dispersion. The effect of matrix diffusion was calculated using an analytical solution to the advection-matrix diffusion problem in which surface retardation was taken into account. Radionuclide transport behavior in rock fractures was explained on the basis of rock structure.     

Mass removal of chlorinated ethenes from rough-walled fractures using permanganate

In situ chemical oxidation (ISCO) employing permanganate is an emerging technology that has been successful at enhancing mass removal from DNAPL source zones in unconsolidated media at the pilot-scale. The focus of this study was to evaluate the applicability of flushing a permanganate solution across two single vertical fractures in a laboratory environment to remove free phase DNAPL. The fracture experiments were designed to represent a portion of a larger fractured aquifer system impacted by a near-surface DNAPL spill over a shallow fractured rock aquifer. Each fracture was characterized by hydraulic and tracer tests, and the aperture field for one of the fractures was mapped using a co-ordinate measurement machine. Following DNAPL emplacement, a series of water and permanganate flushes were performed. To support observations from the fracture experiments, a set of batch experiments was conducted. The data from both fracture experiments showed that the post-oxidation effluent concentration was not impacted by the oxidant flush; however, changes in the aperture distribution, flow field, and flow rate were observed. These changes resulted in a significant decrease to the mass loading from the fractures, and were attributed to the build-up of oxidation by-products (manganese oxides and carbon dioxide) within the fracture which was corroborated by the batch experiment data and visual examination of the walls of one fracture. These results provide insight into the potential impact that a permanganate solution and oxidation by-products can have on the aperture distribution within a fracture and on DNAPL mass transfer rates. A permanganate flush or injection completed within a fractured rock aquifer may lead to the development of an insoluble product adjacent to the DNAPL which results in the reduction or complete elimination of advective regions near the DNAPL and reduces mass transfer rates. This outcome would have significant implications on the plume generating potential of the remaining DNAPL.     

The effects of matrix diffusion on radionuclide migration in rock column experiments

Rock column experiments were performed to examine the effects of matrix diffusion and hydrodynamic dispersion on the migration of radionuclides at the laboratory scale. Tritiated water and chloride transportation was studied in intact mica gneiss and in altered more porous tonalite columns with narrow flow channels. The column diffusion properties were estimated prior to water flow experiments using the gas diffusion method with helium as the tracer gas. The numerical compartment model for advection and dispersion, with and without matrix diffusion, was used to interpret the tracer transport in the columns. Matrix diffusion behavior was also distinguished from dominating hydrodynamic dispersion in rock column experiments at the slowest water flow rates.     

Modeling field-scale multiple tracer injection at a low-level waste disposal site in fractured rocks: effect of multiscale heterogeneity and source term uncertainty on conceptual understanding of mass transfer processes

Multiple factors may affect the scale-up of laboratory multi-tracer injection into structured porous media to the field. Under transient flow conditions and with multiscale heterogeneities in the field, previous attempts to scale-up laboratory experiments have not answered definitely the questions about the governing mechanisms and the spatial extent of the influence of small-scale mass transfer processes such as matrix diffusion. The objective of this research is to investigate the effects of multiscale heterogeneity, mechanistic and site model conceptualization, and source term density effect on elucidating and interpreting tracer movement in the field. Tracer release and monitoring information previously obtained in a field campaign of multiple, conservative tracer injection under natural hydraulic gradients at a low-level waste disposal site in eastern Tennessee, United States, is used for the research. A suite of two-pore-domain, or fracture-matrix, groundwater flow and transport models are calibrated and used to conduct model parameter and prediction uncertainty analyses. These efforts are facilitated by a novel nested Latin-hypercube sampling technique. Our results verify, at field scale, a multiple-pore-domain, multiscale mechanistic conceptual model that was used previously to interpret only laboratory observations. The results also suggest that, integrated over the entire field site, mass flux rates attributable to small-scale mass transfer are comparable to that of field-scale solute transport. The uncertainty analyses show that fracture spacing is the most important model parameter and model prediction uncertainty is relatively higher at the interface between the preferred flow path and its parent bedrock. The comparisons of site conceptual models indicate that the effect of matrix diffusion may be confined to the immediate neighborhood of the preferential flow path. Finally, because the relatively large amount of tracer needed for field studies, it is likely that source term density effect may exaggerate or obscure the effect of matrix diffusion on the movement of tracers from the preferred flow path into the bedrock.     

Field investigation into unsaturated flow and transport in a fault: model analyses

Results of a fault test performed in the unsaturated zone of Yucca Mountain, Nevada, were analyzed using a three-dimensional numerical model. The fault was explicitly represented as a discrete feature and the surrounding rock was treated as a dual-continuum (fracture-matrix) system. Model calibration against seepage and water-travel-velocity data suggests that lithophysal cavities connected to fractures can considerably enhance the effective fracture porosity and therefore retard water flow in fractures. Comparisons between simulation results and tracer concentration data also indicate that matrix diffusion is an important mechanism for solute transport in unsaturated fractured rock. We found that an increased fault-matrix and fracture-matrix interface areas were needed to match the observed tracer data, which is consistent with previous studies. The study results suggest that the current site-scale model for the unsaturated zone of Yucca Mountain may underestimate radionuclide transport time within the unsaturated zone, because an increased fracture-matrix interface area and the increased effective fracture porosity arising from lithophysal cavities are not considered in the current site-scale model.     

Measurements of groundwater velocity in discrete rock fractures

Estimating groundwater velocity in fracture networks using a Darcy or cubic law calculation is complicated by the wide distribution of fracture aperture often found in these systems and by the difficulty in measuring hydraulic head in discrete fracture features. Although difficult to conduct in a fractured rock setting, the point dilution method can be utilized to collect direct measurements of groundwater velocity in individual fractures. To compare measured against calculated velocities, more than 100 point dilution experiments were conducted within a 35 x 35 m area of a single fracture and in discrete fracture features within a fracture network at a larger scale. The dilution experiments were conducted by isolating a fracture feature in a borehole, measuring the hydraulic aperture, and measuring the decay of an injected tracer due to the advective groundwater flux across the fracture. Groundwater velocity was estimated using the hydraulic aperture and the rate of decay of the injected tracer. Estimates of the local hydraulic gradient were calculated via the cubic law using the velocity estimate and the hydraulic aperture. The results of the tests conducted in the single fracture show variable (1 to 33 m/day) but on average higher velocities in comparison to that measured during a natural gradient tracer experiment conducted previously (in which the effects of matrix diffusion were accounted for) and to that which would be calculated using the cubic law. Based on these results, it was determined that the best estimate of the average groundwater velocity, at the scale of the measurement area used for the cubic law calculations, could only be obtained using the largest apertures in the aperture distribution. Variability of the velocity measurements was also observed over time. Increases in velocity were attributed to the effect of rainfall although concurrent increases in hydraulic gradient were not detected (likely within the tolerance of the measuring devices). The groundwater velocities measured in the fracture network varied over a wider range than at the scale of the single fracture (from 2 to 388 m/day). No correlation, however, was observed between the size of the fracture aperture and measured velocity.     

Rigorous uncertainty assessment in contaminant transport inverse modelling: a case study of fluoride diffusion through clay liners

Inverse methods used in assessing landfill liner design have not yet taken advantage of current developments in inverse procedures. Here, a method for inverting contaminant transport models is presented including a general error model and procedures for differentially weighted multiple response regression. General error models are employed in cases where the residuals are heteroscedastic and correlated, and lead to valid inference on model parameter and predictive uncertainty. The Shuffled Complex Evolution algorithm is used to optimise model parameters. Model parameter uncertainty is assessed by exploring the posterior probability distribution with the Metropolis algorithm, a Markov chain Monte Carlo sampling method. The inverse method is applied to simultaneously determine the sorption and diffusion parameters from laboratory diffusion cell experiments. In these experiments, fluoride migration through kaolin clays was measured by sampling the source and collector cells over time. To uniquely determine the transport model parameters, it was necessary to simultaneously fit the observed data from two independent diffusion cell experiments with different initial concentrations. The jointly fitted transport model parameters compared well with those fitted to independent batch experiments.     

Evaluation of tritiated water diffusion through the Toarcian clayey formation of the Tournemire experimental site (France)

Through-diffusion experiments with tritiated water were performed on argillaceous samples from various zones of the Tournemire test site. It was intended to evaluate the homogeneity of the transport property of unfracturated samples and the influence of the orientation and the nature of the samples (presence of an opened fracture or a pre-existing tectonic fracture filled with calcite and pyrite). Homogeneous values of the tritiated water (HTO) effective diffusion coefficients were deduced from experiments carried out when diffusion occurred parallel to the stratigraphic bedding, with an apparent sensitivity to experimental conditions. Anisotropy was significant, De(HTO) perpendicular to the bedding being 1/3 lower than that parallel to the bedding. The observed fractures of the samples created by mechanical stress and partial dehydration during sawing and the presence of a pre-existing opened fracture did not affect the effective diffusion coefficients of tritiated water, which is probably due to the healing ability of the clayey medium during the re-saturation phases of the equilibrium steps performed prior to the diffusion experiments. On the contrary, a significant decrease of this transport parameter was induced by the occurrence of a pre-existing tectonic fracture, which was assigned to the dense structure of the filling phases.     

Matrix diffusion coefficients in volcanic rocks at the Nevada test site: influence of matrix porosity, matrix permeability, and fracture coating minerals

Diffusion cell experiments were conducted to measure nonsorbing solute matrix diffusion coefficients in forty-seven different volcanic rock matrix samples from eight different locations (with multiple depth intervals represented at several locations) at the Nevada Test Site. The solutes used in the experiments included bromide, iodide, pentafluorobenzoate (PFBA), and tritiated water ((3)HHO). The porosity and saturated permeability of most of the diffusion cell samples were measured to evaluate the correlation of these two variables with tracer matrix diffusion coefficients divided by the free-water diffusion coefficient (D(m)/D*). To investigate the influence of fracture coating minerals on matrix diffusion, ten of the diffusion cells represented paired samples from the same depth interval in which one sample contained a fracture surface with mineral coatings and the other sample consisted of only pure matrix. The log of (D(m)/D*) was found to be positively correlated with both the matrix porosity and the log of matrix permeability. A multiple linear regression analysis indicated that both parameters contributed significantly to the regression at the 95% confidence level. However, the log of the matrix diffusion coefficient was more highly-correlated with the log of matrix permeability than with matrix porosity, which suggests that matrix diffusion coefficients, like matrix permeabilities, have a greater dependence on the interconnectedness of matrix porosity than on the matrix porosity itself. The regression equation for the volcanic rocks was found to provide satisfactory predictions of log(D(m)/D*) for other types of rocks with similar ranges of matrix porosity and permeability as the volcanic rocks, but it did a poorer job predicting log(D(m)/D*) for rocks with lower porosities and/or permeabilities. The presence of mineral coatings on fracture walls did not appear to have a significant effect on matrix diffusion in the ten paired diffusion cell experiments.     

Solute transport in crystalline rocks at Aspö--II: blind predictions,inverse modelling and lessons learnt from test STT1

Based on the results from detailed structural and petrological characterisation and on up-scaled laboratory values for sorption and diffusion, blind predictions were made for the STT1 dipole tracer test performed in the Swedish Asp? Hard Rock Laboratory. The tracers used were nonsorbing, such as uranine and tritiated water, weakly sorbing 22Na(+), 85Sr(2+), 47Ca(2+)and more strongly sorbing 86Rb(+), 133Ba(2+), 137Cs(+). Our model consists of two parts: (1) a flow part based on a 2D-streamtube formalism accounting for the natural background flow field and with an underlying homogeneous and isotropic transmissivity field and (2) a transport part in terms of the dual porosity medium approach which is linked to the flow part by the flow porosity. The calibration of the model was done using the data from one single uranine breakthrough (PDT3). The study clearly showed that matrix diffusion into a highly porous material, fault gouge, had to be included in our model evidenced by the characteristic shape of the breakthrough curve and in line with geological observations. After the disclosure of the measurements, it turned out that, in spite of the simplicity of our model, the prediction for the nonsorbing and weakly sorbing tracers was fairly good. The blind prediction for the more strongly sorbing tracers was in general less accurate. The reason for the good predictions is deemed to be the result of the choice of a model structure strongly based on geological observation. The breakthrough curves were inversely modelled to determine in situ values for the transport parameters and to draw consequences on the model structure applied. For good fits, only one additional fracture family in contact with cataclasite had to be taken into account, but no new transport mechanisms had to be invoked. The in situ values for the effective diffusion coefficient for fault gouge are a factor of 2-15 larger than the laboratory data. For cataclasite, both data sets have values comparable to laboratory data. The extracted K(d) values for the weakly sorbing tracers are larger than Swedish laboratory data by a factor of 25-60, but agree within a factor of 3-5 for the more strongly sorbing nuclides. The reason for the inconsistency concerning K(d)s is the use of fresh granite in the laboratory studies, whereas tracers in the field experiments interact only with fracture fault gouge and to a lesser extent with cataclasite both being mineralogically very different (e.g. clay-bearing) from the intact wall rock.     

Characterization of fracture aperture field heterogeneity by electrical resistance measurement

We use electrical resistance measurements to characterize the aperture field in a rough fracture. This is done by performing displacement experiments using two miscible fluids of different electrical resistivity and monitoring the time variation of the overall fracture resistance. Two fractures have been used: their complementary rough walls are identical but have different relative shear displacements which create "channel" or "barrier" structures in the aperture field, respectively parallel or perpendicular to the mean flow velocity U(→). In the "channel" geometry, the resistance displays an initial linear variation followed by a tail part which reflects the velocity contrast between slow and fast flow channels. In the "barrier" geometry, a change in the slope between two linear zones suggests the existence of domains of different characteristic aperture along the fracture. These variations are well reproduced analytically and numerically using simple flow models. For each geometry, we present then a data inversion procedure that allows one to extract the key features of the heterogeneity from the resistance measurement.     

Flamelet modeling of NO formation in laminar and turbulent diffusion flames

The applicability of the laminar flamelet concept for the formation and destruction of nitric oxides in laminar and turbulent diffusion flames has been studied. In a first step, temperatures and species concentrations in an axisymmetric laminar diffusion flame have been calculated (i) by solving the detailed conservation equations and (ii) by applying the laminar flamelet concept. The main purpose of this step was the identification of differences between results from both approaches. It turned out that for highly temperature sensitive or relatively slow chemical processes, the inclusion of the full range of the prevailing scalar dissipation rates plays a major role for the calculated species concentrations. This behavior is obvious from the concept of the laminar flamelet model, where the scalar dissipation rate can be discussed in terms of the reciprocal of a residence time for attaining chemical equilibrium. In a second step, flamelet modeling of NOx formation was extended to a turbulent hydrogen diffusion flame. In both the steps, the flow fields of the flames were calculated by solving the Navier-Stokes equations in axisymmetric formulation using the SIMPLER algorithm. For the turbulent flow, Favre-averaged equations have been used and turbulence was modeled with the standard k-epsilon model including a correction term for axisymmetric systems. The averaging of the species concentrations was accomplished with presumed shape probability density functions (pdfs). The pdf of the mixture fraction was described with a beta-function whereas that of the scalar dissipation rate was assumed to be log-normal. Buoyancy effects have been taken into account. The calculated temperatures and concentrations were compared with data from different experiments.     

Predictive modelling of NAPL injection tests in variable aperture spatially correlated fractures.

In preparation for a field experiment where a NAPL will be injected into a fractured sandstone aquifer, a 2D invasion percolation model of DNAPL migration in a single horizontal fracture with varying aperture has been developed. This simulation investigated the effect of spatially correlated fracture aperture fields on pressure-saturation relationships as a function of variable aperture mean, standard deviation, and spatial correlation statistics under hydrostatic conditions. Results from spatially correlated variable aperture fields can be significantly different from random fields. Longer ranges decreased entry pressures and higher standard deviations decreased nonwetting phase saturations. Mean aperture is the major control on displacement pressure, entry pressure and the form of the pressure-saturation curve. Simulation results using statistical parameters for a variable aperture natural sandstone fracture as described by Yeo et al. [International Journal of Rock Mechanics and Mining Sciences 35 (1998) 1051] closely resemble a Brooks-Corey pressure-saturation function, and exhibit a distinct entry pressure followed by a rapid increase in nonwetting phase saturation. Graphical estimates of entry pressure provide a good approximation of the critical aperture controlling the formation of a continuous nonwetting phase pathway in a variable aperture fracture. Consequently, we show that multiphase flow concepts developed for porous media can successfully be applied to variable aperture fractures. Entry pressure correlates well to the mean aperture in these simulations when using a Gaussian distribution of fracture aperture. Interpreted aperture distributions from NAPL injection experiments do not fit the true distribution well at low nonwetting phase saturations, but do at higher saturations above the entry pressure. Consequently, true, mechanical aperture variation within a fracture plane cannot be determined from NAPL injection experiments either in the field or laboratory.