## Abstract

This paper presents a technique for determining the near surface permeability of geomaterials and involves the application of a uniform flow rate to an open central region of a sealed annular patch on an otherwise unsealed flat surface. Darcy’s flow is established during attainment of a steady pressure at a constant flow rate. This paper describes the experimental configuration and its theoretical analysis via mathematical and computational techniques. The methods are applied to investigate the surface permeability characteristics of a cuboidal block of Indiana limestone measuring 508 mm. An inverse analysis procedure is used to estimate the permeability characteristics at the interior of the Indiana limestone block. The resulting spatial distribution of permeability is used to estimate the *effective permeability* of the tested block.

## 1. Introduction

Permeability is a key parameter in environmental geosciences and geomechanics problems dealing with groundwater hydrology (Bear 1972; de Marsily 1986), groundwater contamination (Bear *et al.* 1993; Selvadurai 2003*a*, 2004a, 2008), geothermal energy extraction (Dickson & Fanelli 1995; Duffield & Sass 2003), geological sequestration of CO_{2} (Bachu & Adams 2003; Moore *et al.* 2004; Tsang *et al.* 2008), deep geological disposal of contaminants (Apps & Tsang 1996; Selvadurai 2006) and geological storage of heat-emitting nuclear fuel wastes (Laughton *et al.* 1986; Selvadurai & Nguyen 1997). The property of permeability or intrinsic permeability depends solely on the accessible pore structure to fluid flow in a geomaterial (Brace 1977; Bernabé *et al.* 1982). Particularly in hydrogeology, permeability is assumed to be an unchanging property, whereas geomechanical factors including the stress state can result in the generation of defects in the intact geomaterial, which can significantly influence its permeability (Heidland 2003; Selvadurai 2004*b*; Selvadurai & Głowacki 2008). Reliable estimation of permeability of geomaterials is crucial to computational models that couple flow, deformation and transport processes (Bear *et al*. 1993; Selvadurai & Nguyen 1995; Stephansson *et al.* 1996; Alonso *et al.* 2005).

While scaling approaches have been proposed in the literature (Neuman & di Federico 2003), there are no universally accepted guidelines that can be used to identify the relevant scales of interest to permeability estimation; these can range from the crustal level involving dimensions of the order of 0.5–5.0 km, to borehole scales of the order of 30–300 m, to laboratory scales of the order of 5–15 cm. The estimation of ‘permeability’ of a geomaterial can also be significantly influenced not only by natural inhomogeneities that include fractures, fissures, damage zones, voids and vugs but also by coupled geomechanical, geochemical and geomorphological processes. Conventional laboratory measurement of permeability invariably relies on cylindrical samples that have substantially smaller dimensions (maximum 10 cm in diameter and 20 cm in length) in comparison with the dimensions of the geomaterial volume it is expected to represent in a subsurface fluid flow analysis problem. The use of larger specimens can eliminate some of the drawbacks of small-scale testing but the use of conventional axial or radial flow tests with large samples is non-routine. This paper discusses the possible use of a mini-permeameter to determine the local-scale permeability. Ideally, such a device should be versatile enough to allow estimation of local-scale permeability both rapidly and accurately. When the large sample being tested has a plane surface, it is possible to use a patch permeability test, where an annular region provides an impervious seal and the central cavity forms the pressurized zone that establishes steady flow. Although these non-invasive tests can be time-consuming, they present a viable approach in situations where the overall permeability of the block sample is needed in a non-destructive manner. This paper focuses on the development of an experimental procedure that can be used to estimate the surface permeability distribution in a cuboidal sample of Indiana limestone measuring 508 mm. The term ‘surface permeability’ is intended to signify the permeability of an effective support volume of the porous medium, which is determined by initiating steady flow at an aperture located at a partially sealed surface. In this research the dimensions of the support volume depend on the external diameter of the annular sealing region.

Laboratory procedures developed for steady-state permeability testing of rock cores invariably involve the application of steady-state flow, usually in the longitudinal direction of the cylindrical sample. The main advantage of a steady-state test is that the interpretation of the results is relatively straightforward, dependent only on the hydraulic boundary conditions of the test, the dimensions of the sample and a knowledge of the steady flow rate. In materials with relatively low permeability, an inordinate amount of time is required to attain steady-state flow conditions, and recourse is usually made to transient methods such as hydraulic pulse tests. Research dealing with the evaluation of permeability properties of geomaterials using both steady-state and pulse tests is quite extensive; references to seminal work and current literature can be found in the article by Selvadurai (2009) (see also Brace *et al.* 1968; Hsieh *et al.* 1981; Selvadurai & Carnaffan 1997; Selvadurai *et al.* 2005; Selvadurai & Selvadurai 2007; Selvadurai & Głowacki 2008).

Large samples, a prerequisite to conducting surface permeability tests, can usually be obtained from either rock outcrops or from tunnels adits and test pits that are used for other geological investigations. The test involves the application of a constant flow rate to the central opening of an annular sealed section of the surface of the test specimen, such that steady flow conditions are attained. The use of localized steady flows for permeability testing has been reported by Dykstra & Parsons (1950) and these experimental techniques were rigorously modelled by Goggin *et al.* (1988). Innovative research along these lines was also conducted by Tidwell & Wilson (1997) in connection with the measurement of air permeability of a large block of Berea sandstone and the errors involved in the numerical evaluation of the intake shape factors were examined by Tartakovsky *et al.* (2000). This paper briefly discusses the experimental procedure, provides a mathematical treatment of special cases of the experimental configuration and uses computational models to interpret the data for general situations. Finally, representative results for the surface permeability determined by conducting patch permeability measurements at nine locations on each surface of the cubical block of Indiana limestone are given. The experimental results, along with estimates for the interior permeability distribution derived from kriging procedures and theoretical relationships proposed in the literature, are used to determine the *effective permeability* of the entire block.

## 2. The test facility and experimental procedures

Indiana limestone has been used quite extensively in research investigations, and the geomechanical (skeletal Young’s modulus and Poisson’s ratio; tensile and compressive strengths), physical (porosity, unit weight) and mineralogical characteristics of the rock are well documented in the open literature (Głowacki 2007). The faces of the sample were saw-cut, creating a surface texture consistent with that of an FEPA grade P120 emery paper. Actual surface roughness characterizations were not required for the research investigation.

The test required the development of a perfect hydraulic seal over an annular region in contact with the rough surface of the cuboidal specimen and the application of a constant flow rate to the exposed central region, to attain steady-state flow. The dimensions of the annular sealing region can vary depending on the experimental configuration, but the attainment of a perfect seal over the annular region is essential for the success of the experiments. Figure 1 shows a cross section through the permeameter and a schematic view of the experimental facility is shown in figure 2. A reaction frame was used to provide the sealing loads and a manually operated hydraulic cylinder and a load cell (interface 10 000 lb (approx. 4545 kg) loadcell; accuracy ±0.05%) was used to maintain the sealing loads. The LC8A Shimadzu pump provided a prescribed flow rate accurate to within ±2 per cent of the prescribed flow rate, with a maximum pressure range of 0–30 MPa. The accuracy of the pump was also verified separately through direct measurement of the flow rate. The permeameter is fitted with a pressure transducer (DCT Instruments PZA 100 AC; 0–689.4 kPa with an accuracy of ±0.05%) and a port for de-airing the cavity region. All instrumentation was monitored using a Techmatron Instruments Inc. Data Acquisition System. Water used was regular tap water, which was not de-aired or de-mineralized but filtered through a 10 μm filter and maintained at a constant room temperature of 23^{°}C. The limestone block was kept immersed in a water reservoir maintained at room temperature. The cross-head of the reaction frame could be moved to accommodate specific test locations on the surface of the sample, as shown in figure 3. This allowed measurement of the distribution of surface permeability on all faces of the cuboidal Indiana limestone specimen.

### (a) The permeameter

Sealing between the permeameter (figure 1) and the rock surface was achieved by applying a load to a natural rubber gasket (40 Shore A durometer) with an inside diameter of 25.4 mm and outside diameter of 101.6 mm. The interior to exterior radii (or diameter) ratio (*a*/*b*) is also referred to as the ‘tip seal ratio’ and in these experiments it was set at 4. The sealing annulus was confined internally and externally by retaining rings to prevent in-plane radial expansion, both during the application of the sealing load and during pressurization of the interior cavity. The sealing technique is efficient and less elaborate than that used by Tidwell & Wilson (1997) and is described elsewhere (Selvadurai, P. A. 2010). The sealed cavity was de-aired by flushing the system with water.

### (b) Sealing procedure

In these experiments, we adopted the procedure suggested by Tidwell & Wilson (1997) for determining an adequate sealing load. A separate experiment was conducted where the sealing load was increased in stages and the permeability measured at the same location. The observed reduction in permeability is directly related to the compression of the gasket to conform to the surface topography of the limestone, which prevents fluid migration through the interface. At gasket compression stresses above 1.4 MPa, there were no appreciable changes (i.e. less than 2%) in the estimated permeability. The mechanical effects of pore structure alteration owing to sealing pressures of up to 2.5 MPa are considered to have a negligible influence, and studies by Selvadurai & Głowacki (2008) indicate that confining pressures greater than 10 MPa are needed to induce appreciable alterations in permeability owing to pore closure, pore collapse and other irreversible pore structure alterations. To account for any non-uniformity of the surface texture at the contact locations, the sealing pressure was maintained at 1.75 MPa for all test locations. The effectiveness of the sealing stress was also established by repeating the test at the end of the series of experiments.

### (c) Test procedure and results

If the permeability is to be defined by an appeal to Darcy’s Law, then the induced flow should satisfy conditions applicable to the Darcy flow regime. Extensive discussion of the limits of applicability of Darcy flow in porous media and references to further studies are given by Zeng & Grigg (2006). The two criteria—the Reynolds number and the Forchheimer number—are in general used for assessing limits to the applicability of Darcy flow. The Reynolds number *Re*=*ρ*_{f}*vd*/*μ*, where *ρ*_{f} is the fluid density, *v* is the superficial flow velocity, *d* is a characteristic length defined in relation to either the mean pore throat size or the mean grain size (Hornberger *et al.* 1998) and *μ* is the dynamic viscosity. The Forchheimer number *Fo*=*kβρ*_{f}*v*/*μ*, where *k* is the permeability at ‘zero’ velocity and *β* is a non-Darcy coefficient. The latter criterion is less than helpful if the research focuses on the determination of the permeability of the porous medium. A simpler estimate for limits of applicability of Darcy flow particularly suitable for porous rocks is given by Philips (1991) and it requires *Re*≪*n*, where *n* is the porosity. The measurement of the pore throat diameter in porous rocks is a non-routine exercise. The experimental results given by Wardlaw *et al.* (1987) indicate that the pore throat diameter for Indiana limestone can vary from a minimum of 2×10^{−9} mm to a maximum of 46×10^{−9} mm, with an average of 11×10^{−9} mm. Churcher *et al.* (1991) indicated a pore throat diameter of approximately 15×10^{−5} mm. The porosity of Indiana limestone measured from water saturation tests (Selvadurai & Głowacki 2008) was approximately 0.17. Considering an average superficial flow velocity of *v*≈3.3×10^{−4} m s^{−1} (consistent with a peak flow rate of 10 ml min^{−1} over the central opening of 25.4 mm), and on the basis of these pore throat diameter estimates, we obtain *Re*≈(3×10^{−9},5×10^{−5}), which is significantly smaller than the estimate for porosity. On the basis of this estimate, we conclude that all the experiments performed were within the Darcy flow range. The validity of the Darcy flow regime was further confirmed by performing variable flow rate tests and observing the deviations from linearity between flow rate and steady-state pressure (figure 4). The peak flow rate of 10 ml min^{−1} is well within the flow regime required for the applicability of Darcy flow. In this paper, results of 162 steady-state permeability tests conducted on the six faces of the block of Indiana limestone are reported; a minimum of three tests were conducted at a known temperature at each of the nine locations per side (figure 3), and the tests were interpreted using mathematical results and computational approaches described in the ensuing sections. At each location, the maximum and minimum pressures reached during attainment of a steady-state flow were recorded. Previous experiments (Selvadurai & Głowacki 2008) conducted on cylindrical cores of Indiana limestone (100 mm diameter and 200 mm long) subjected to one-dimensional flow indicate little or no mineral dissolution that could influence permeability alterations. Similarly, the surface permeability tests also indicate virtually no alteration in the estimated value of the permeability in subsequent tests. The maximum variability of approximately 2.28 per cent during tests conducted at the same location within a five-month interval is within the range of sensitivity of the experimentation and cannot be attributed to any physical alterations of the pore structure.

## 3. Mathematical modelling of the permeameter

To estimate the permeability from the results of steady-state flow tests, we require appropriate mathematical and/or computational relationships between the measured responses (cavity pressure and flow rate) and the geometry of the annular region. In its most general form, Darcy’s Law can be written as (Bear 1972; Selvadurai 2000*a*)
3.1
where ** v**(

**) is the velocity vector,**

*x***(**

*K***) is a permeability tensor, which is a symmetric, positive definite, second-order tensor,**

*x**Φ*(

**) is the reduced Bernoulli potential consisting of the pressure and datum heads,**

*x**μ*is the dynamic viscosity of water,

*γ*

_{w}is its unit weight and ∇ is the gradient operator. While the theoretical concepts in this definition are relatively straightforward, the determination of the permeability characteristic for such media through the consideration of generalized inverse problems is a more complicated exercise in both mathematical and experimental modelling. A formulation of the non-local aspects of permeability and its relationship to the statistical heterogeneity of the pore structure have been considered by several investigators (Renard & de Marsily 1997; Tartakovsky 2000; Zhang 2002; Neuman & di Federico 2003). The approach adopted here is to assume that the surface permeability test will allow the estimation of the equivalent isotropic permeability of the representative volume or support volume that is addressed by the permeability test and to use the support volume-wise surface permeability values to develop a pattern of permeability distribution for the entire cuboidal limestone specimen, including its interior. Even with this simplification, the estimation of the permeability characteristics of the interior of the cuboidal block requires further analysis. With this limitation in mind, we proceed to make the assumption that the representative volume tested with the surface permeability test can be approximated by an equivalent isotropic and homogeneous form of Darcy’s Law 3.2 where

*K*is the scalar permeability. With this form of Darcy’s Law, for incompressible flow in the pore space and for non-deformability of the porous medium, the steady-state flow can be described by Laplace’s equation 3.3

The sealed annular patch in the surface permeability test can be located at any position on the face of the cuboidal limestone block, making the problem solvable only by appeal to computational solutions of equation (3.3). There are, however, situations that lend themselves to analytical treatment, provided certain assumptions can be made with regard to the dimensions of the annular sealing patch relative to the overall flow domain. The first of these relates to the case of the annular sealing patch located at the centre of the block. In this case, the problem is nearly *axisymmetric* and can be conveniently formulated, analytically, if the porous domain is also assumed to be semi-infinite, i.e. equation (3.3) takes the form
3.4

A Hankel transform development (Sneddon 1966; Selvadurai 2000*b*) gives the solution
3.5
applicable to *z*≥0, where *φ*(*ξ*) is an arbitrary function and *J*_{0}(*ξr*) is the zeroth-order Bessel function of the first kind. The form of the function (3.5) ensures that as . The arbitrary function *φ*(*ξ*) needs to be determined by solving a mixed-boundary value problem applicable to the permeability test. Considering the annular sealed patch problem, we assume that at steady-state: (i) a constant potential *Φ*_{0} is applied to the inner region (0≤*r*≤*a*), (ii) the sealing annular patch (*a*<*r*<*b*) is subjected to null Neumann boundary conditions, and (iii) the region exterior to the annular patch is subjected to zero potential. These boundary conditions yield a set of triple integral equations of the form
3.6
3.7and
3.8
These triple integral equations have been studied in connection with potential theory (Tranter 1960; Cooke 1963; Williams 1963) and classical elasticity problems dealing with annular punch, crack and inclusion problems (Gubenko & Mossakovski 1960; Gladwell 1980; Selvadurai & Singh 1984*a*,*b*). Approximate solutions to the triple system were presented by Goggin *et al.* (1988) and, in a subsequent analysis, Tartakovsky *et al.* (2000) provided an informative discussion of the errors and accuracy involved in the development of numerical solutions. The objective here is to develop an analytical solution, based on the series approximation techniques, to the triple integral equations along the lines similar to those proposed by Selvadurai & Singh (1984*a*,*b*) and Selvadurai (1996) and summarize the results for completeness. We assume that *φ*(*ξ*) admits a representation
3.9
and these unknown functions can also be represented in the forms
3.10and
3.11
These representations give rise to a set of coupled Fredholm integral equations of the form
3.12and
3.13
where *s*=*aζ* and *c*=*a*/*b*. The set of coupled integral equations (3.12) and (3.13) can also be solved using a variety of numerical schemes (Atkinson 1997); here, we assume that the functions *F*_{1}(*ζ*) and *F*_{3}(*ζ*) can be expanded as power series in terms of a small non-dimensional parameter, which is chosen as *c*<1, i.e.
3.14
Expanding (*η*^{2}−*c*^{2}*ζ*^{2})^{−1}and (*ζ*^{2}−*c*^{2}*η*^{2})^{−1} also in power series in terms of the parameter *c* and substituting these into the coupled integral equations (3.12) and (3.13), we can determine the functions and . The potential flow problem for the annular sealed region of the halfspace is thus formally solved, albeit in an approximate series form.

When the outer boundary of the annular region , the results (3.6)–(3.8) yield a classical two-part mixed boundary value problem in potential theory. The result for the potential problem where the circular cavity is subjected to a potential *Φ*_{0} and the region exterior to it is impervious is given by
3.15
The result of particular interest to the analysis of the experiment is the relationship between the potential *Φ*_{0} and the steady flow rate, *Q*, i.e.
3.16
where *S*_{D} is the entry region. Using equations (3.2) and (3.15) in equation (3.16), we obtain
3.17
where *F*=4. The result (3.17) is presented in a form similar to that proposed in geomechanics (Selvadurai 2003*b*) to identify the *intake shape factor* for the entry point, which depends on the geometric characteristics of the intake. In equation (3.17), however, *F* is non-dimensional.

The second limiting case deals with the situation where the region 0≤*r*≤*a* is subjected to the potential *Φ*_{0} and the exterior region is subjected to zero potential. In this case, the potential *Φ*(*r*,*z*) is given by
3.18
We can calculate the flow rate through the pressurized circular patch by using equations (3.18), (3.2) and (3.16); this gives the flow velocity distribution on the region as
3.19and
3.20
Evaluating the velocity components on the plane *z*=0, we note that *v*_{r}(*r*,0)≡0 for 0≤*r*<*a*, and the axial velocity in the circular region is given by
3.21
which is an integral of the Lipschitz–Hankel type, extensively investigated by Eason *et al.* (1955), who also provided numerical results for a large class of such integrals and discussed their convergence. Although the integral equation (3.21) can be expressed in terms of elliptic integrals, the steady-state flow as defined by equation (3.16) is divergent since the function has a non-integrable singularity at *r*=*a* (e.g. Salamon & Walter 1979). The approximate solution for the system of triple integral equations yields, for *c*=(*a*/*b*)≪1, the following result for the flow rate:
3.22
where
3.23
The expansion equation (3.23) will be applicable only for a limited range of the values of *c*. Figure 5 illustrates the variation in the intake shape factor *F*(*c*) as a function of (*b*/*a*) and the excellent agreement with results given by Goggin *et al.* (1988) and Tartakovsky *et al.* (2000) is noted (accurate for values of *a/b* of approx. 0.833).

The extension of the problem to anisotropy of the porous medium entails a prohibitive amount of analysis, which will not be conducive in estimating, in a unique fashion, the full permeability tensor at the location where the test is carried out. A simpler approach is to assume transverse isotropy of the porous medium with coordinate directions aligned with the principal directions (e.g. Tartakovsky *et al.* 2000; Selvadurai 2003*b*, 2004*c*). This analysis assumes that the plane of transverse isotropy is parallel to the boundary plane of the halfspace, which again constrains the interpretation of the data. Even with assumptions of transverse isotropy, the unknowns in the study of the porous medium flow problem include the two principal values of permeability and the inclination of the plane of transverse isotropy to the boundary plane of the halfspace region, and three independent tests are needed to estimate the unknowns, uniquely. Analytical solutions that can examine the problem of the annular sealed area of the surface of a halfspace region are not available. Only recently has a solution been developed (Selvadurai in press) for the case where the circular intake is located at the boundary of a porous medium with transverse isotropy and accounts for inclination (*α*) of the planes of stratification to the boundary plane. The relevant result is given by
3.24
where , *K*_{n} and *K*_{t} are, respectively, the permeabilities normal and along the principal planes of transverse isotropy, and
3.25
In equation (3.25), *K*(*ρ*) is the complete elliptic integral of the first kind, defined by
3.26

As is evident, the circular entry point can provide only limited information that is insufficient to determine the complete structure of the permeability tensor for a transversely isotropic porous medium. The entry point geometry needs to be altered to provide independent interpretations of the influence of transverse isotropy on the flow rates.

## 4. Computational modelling

When considering the cuboidal Indiana limestone specimen used in the experiments, the analytical solution equation (3.22) for the flow rate at the annular sealing patch is most suitable for the central test locations. However, the analytical results are essential for the validation of computational approaches, which will be used to estimate permeability at other locations. It should also be noted that results obtained by Goggin *et al.* (1988) and Tartakovsky *et al.* (2000) indicate that the halfspace solution is even applicable to cylindrical cores that are partially sealed by an annular patch on a plane surface and subjected to null Dirichlet boundary conditions on the remaining surfaces, which represents an exterior domain of infinite permeability. Several computational approaches including boundary element techniques (Gaul *et al.* 2003) and finite-element approaches (Zienkiewicz & Taylor 2006) can be used to solve Laplace’s equation. The boundary element technique can accurately model singularities in the velocity field that can be encountered at a Dirichlet–Neumann boundary and give robust, convergent solutions for finer mesh division but is less effective in modelling inhomogeneities. The domain techniques are straightforward Galerkin approaches and special elements have to be incorporated to rigorously model singular fields (Aalto 1985) or infinite domains (Selvadurai & Karpurapu 1989). Also, domain methods are capable of solving problems for inhomogeneous regions quite effectively. The computational modelling techniques are essential for determining the flow rates for permeameter positions remote from the central location. The required analysis can be performed quite conveniently using one of the several commercially available finite-element codes; here, we have chosen the Comsol multiphysics code in view of further applications involving coupled thermo-poroelastic responses using the same experimental facilities.

### (a) Validation of the computational modelling

To establish the accuracy of the finite-element approach, we consider the potential problem for a semi-infinite domain, where a constant potential (*Φ*_{0}) Dirichlet boundary condition is prescribed over the region 0≤*r*<*a*, and a null Neumann boundary condition is prescribed over on the plane *z*=0. The results of the finite-element computations were compared in relation to the exact result for the flow rate given by equation (3.17). Two types of finite-element discretizations that employ tetrahedral elements were used to assess the accuracy of the computational approach. In the first, the flow domain is zoned into three regions of varying mesh refinement and 299 945 elements (figure 6*a*), and, in the second, a single zone with a finely graded mesh in the vicinity of the entry point resulting in 252 227 elements was used (figure 6*b*). A robust solver SPOOLES, provided in Comsol, was used to generate the steady-state solution. The computations using these models were able to predict the permeability (for a given flow rate *Q* and flow potential *Φ*_{0}) to within 2.74 per cent of the analytical result using the mesh discretization shown in figure 6*a* and 2.05 per cent using the mesh discretization shown in figure 6*b*. The discrepancy can largely be attributed to the absence of singularity elements but, more importantly, to the absence of infinite elements in the computational approach. The advantages that can be gained from incorporating such elements are not apparent, particularly in view of the relatively small errors involved. We can conclude that, when the outer radius *R*_{D} of the domain being modelled is larger than 8*a* and the length *L*_{D} is larger than 40*a*, the solution for the annular patch seal with *b*/*a*=4 can be approximated by the corresponding result for the circular entry point solution associated with a halfspace region, the maximum error being less than 1.95 per cent. Computations were also performed to compare the results obtained for experimental configurations, where the permeameter was located either along the edge or at the corner of the cuboidal block. The symmetries associated with these configurations can be used to reduce the size of the domain that needs to be modelled. The reduced domains and the relevant finite-element discretizations are shown in figure 7. The computational estimates for the case where the permeameter is located at the edge under-predicts the analytical result given by equation (3.17) by 3.07 per cent and the corresponding result for the permeameter located at the corner of the block under-predicts the analytical solution by 3.44 per cent. In summary, the computational calibrations performed indicate that the analytical solution for the circular entry point located at the surface of an impervious halfspace can be used to model the permeameter with a tip seal radii ratio (*b*/*a*)>4, to within an accuracy of approximately 3.5 per cent. The computational modelling will, however, be used to estimate the permeability in the vicinity of the sealing region of the permeameter. The issue of the ‘blind spot’ of the permeameter, where the domain exterior to a toroidal region (for a homogeneous isotropic porous medium) remains unprobed by the permeameter, has been discussed by Tartakovsky *et al.* (2000) and Molz *et al.* (2003). Since the extent of the sampled porous medium depends on the permeability properties that are being sought, the issue cannot be resolved with certainty; suffice it to note that, in general, for a hydraulically homogeneous and isotropic porous medium, the permeameter will sample a domain roughly equal to a hemispherical region of radius equal to the outer radius *b*. The results obtained by Goggin *et al.* (1988) indicate that, when (*b*/*a*)>4 and if isotropy is assumed, the intake shape factor for the annular region with a centrally pressurized cavity is almost identical to the circular intake located at an impervious halfspace.

### (b) Computational modelling of a subregion

The basic premise of the analysis of the surface permeability test is that (i) it examines the effective local permeability, which can be regarded as an estimate for the equivalent homogeneous permeability associated with the support volume, similar to the hemispherical region described previously, and (ii) it is relatively uninfluenced by permeability values exterior to the support volume. To establish the accuracy of this assertion, computational modelling was performed on a subcube region measuring 127 mm and sealed at the surface by an annular patch of internal radius 12.7 mm and external radius 50.8 mm. Referring to figure 8, the surface *S*_{0} is subjected to a constant potential *Φ*_{0}, the annular sealing *S*_{1} was subjected to null Neumann conditions, the region *S*_{2} was subjected to null Dirichlet conditions and the plane surfaces *S*_{3} and *S*_{4} were subjected to either (i) null Dirichlet boundary conditions or (ii) null Neumann boundary conditions. The relevant boundary value problems were examined using the Comsol computational model, as shown in figure 8. For the case where the surfaces *S*_{3} and *S*_{4} were subjected to null Dirichlet boundary conditions, the ratio between the computational result and the halfspace solution is 1.010. For the case where the surfaces *S*_{3} and *S*_{4} were subjected to null Neumann boundary conditions, the ratio between the computational result and the halfspace solution is 0.9557. Hence, the maximum error associated with even extreme limits of permeability external to the local cuboidal domain is less than 5.43 per cent, which is regarded as an acceptable margin of error in the estimation of permeability.

In addition, an ideal cuboidal region was constructed using a random set of 64 subcubes, where the permeabilities varied between the experimental range of 250×10^{−15} m^{2} and 11×10^{−15} m^{2}. A subcube chosen at a random location yields a local permeability (estimated from equation (3.22)) of 152.62×10^{−15} m^{2}. The same subcube examined in isolation, using the approach described here, gives permeabilities between 148.31×10^{−15} m^{2} and 156.52×10^{−15} m^{2}, corresponding to the two types of boundary conditions. It is clear that the average permeability is a representative measure of the local value of the tested region and largely uninfluenced even by extreme variations in permeabilities beyond the zone of influence of the permeameter.

## 5. Experimental results

Each face of the cuboidal Indiana limestone sample was tested at nine locations and a minimum of three tests were conducted at each location to determine the local steady-state central region pressure consistent with an assigned flow rate. The steady-state pressure stabilized within approximately 30 min of application of the steady-state flow rate (figure 4). The variability in the steady-state pressures for succesive tests conducted at any given location was approximately 3.58 per cent. Table 1 illustrates the average values for the permeability at the 54 test locations as determined from the computational modelling. Figure 9 illustrates the initial data for the 54 locations, where the area over which the permeability is assigned corresponds to the outer radius of the sealing region.

An important aspect of the research is to completely use the results of the surface permeability measurements to determine plausible variations of permeability at the *interior* of the block of Indiana limestone. There are several approaches that can be used to extrapolate the data to generate the interior mappings of permeability variation and, in this investigation, we employ a ‘kriging’ technique. Kriging techniques have been used quite successfully in a wide range of scientific and engineering endeavours, including oceanography, biosciences, neurosciences, nanotechnology, mineral resources engineering and materials science, where measurements derived on the surface of a region are used to map the distributions in the interior of the domain. The mathematical basis for the procedure is well established (Journel & Huijbregts 1978; Deutsch & Journel 1997) and accounts of the application of kriging techniques in hydrogeology are given in the treatise by de Marsily (1986) and in the text by Kitanidis (1997). Furthermore, since the surface permeability estimates are unchanged during the experiments (i.e. no mineral dissolution, stress-induced alterations of the permeability, etc.), the kriging technique has a greater ability to predict the permeability distribution at the interior of the block. For this purpose, the EasyKrig v.3.0 available in the Matlab Kriging Toolbox was used. The method uses variograms to express spatial variations and minimize the error of the predicted values, which are estimated by the spatial distribution of the known values. Figure 10 shows the dataset as imported into Matlab for the kriging application. Kriging can be conducted using various models, of differing levels of accuracy, that are characterized by their ability to correctly simulate the experimental variogram. The five tested models, in increasing order of accuracy, were *linear* (*Q*_{1}=−0.1535, *Q*_{2}=1.602), *spherical* (*Q*_{1}=−0.1393, *Q*_{2}=1.456), *normalized sinc function* (*Q*_{1}=0.0301, *Q*_{2}=0.7873), *Gaussian* (*Q*_{1}=−0.1009, *Q*_{2}=0.8983) and *exponential Bessel* (*Q*_{1}=−0.0919, *Q*_{2}=0.9865). For this reason, the chosen model was an ‘exponential Bessel’ variogram as seen in figure 11. Kitanidis (1997) has shown that this model with the above values of *Q*_{1} and *Q*_{2} satisfies statistics for the sampled data, thus implying that there is only a 5 per cent probability that the correct model has been excluded. The kriging procedure can be regarded as a ‘smoothing’ process, which generally tends to reduce the higher values of permeability and vice versa. EasyKrig 3.0 permits visual observation of the procedure to ensure that changes made to the permeability through the kriging algorithm are within a prescribed tolerance. The average change in permeability from the initial data to those based on the kriging procedure was (0.01049)×10^{−15} m^{2}. The kriging estimates for permeability distributions within the Indiana limestone block, as estimated from surface permeability measurements, are shown in figure 12. The maximum or minimum estimates for the permeabilities indicated in table 1 gave kriged data almost identical to the results shown in figure 12.

## 6. The effective permeability of the cuboid

It is pertinent to examine the results of the research in the context of the scale at which the results can be applied in a practical context; this leads us to the evaluation of the *effective permeability* that can be assigned to the tested cuboidal sample of Indiana limestone. Despite the heterogeneity that is inherent in geological formations of sedimentary origin, the effective permeability is regarded as a quantity that can indicate the effective fluid transport characteristics. The effective permeability of a region can thus be defined as the permeability that gives rise to the same flux or flow rate, under the same hydraulic gradient applied to the region. For example, if the dimensions of the porous medium being examined are on the same scale as the cuboidal region that was tested, then the permeability distributions determined from the point estimates and kriging procedures can be regarded as the actual values that need to be used in flow calculations. When examining the data in the context of a region that is much larger than the dimensions of the block that was tested (say several kilometres), the block then corresponds to a mathematical point, which, in general, could exhibit hydraulic anisotropy (with spatial dependency being excluded by definition). The topic of calculation of effective permeability of a heterogeneous porous medium has been the subject of extensive research, particularly from theoretical perspectives. Comprehensive expositions can be found in the articles by Cushman (1986), Dagan (1989, 1993), Noetinger (1994), Christakos *et al.* (1995), Wen & Gómez-Hernández (1996), Renard & de Marsily (1997) and Markov & Preziosi (2000). These articles and volumes contain, *inter alia*, a large number of references to research on this topic. The earliest of the rigorous treatments involving statistical continuum theories of heterogeneous media with random geometry was given by Beran (1968). Theories have also been developed by a number of authors (e.g. Batchelor 1974) for the estimation of the *effective permeability* of highly porous media. Of related interest is the discussion of the concept of a representative volume element (RVE) pertaining to Darcy flow in a random medium examined by Du & Ostoja-Starzewski (2006). These authors have applied the self-consistent theories developed by Hill (1963) for heterogeneous elastic media and two-dimensional computational approaches to establish criteria for limits of validity of the RVE on the relative pore geometry. These advances are perhaps not directly applicable to a range of porosities encountered in the Indiana limestone that was tested (i.e. there was no compelling visual evidence of dispersion of inhomogeneous permeable inclusions in a permeable matrix with widely different permeabilities (e.g. McCarthy 1991)). The work of King (1987, 1989) indicates that, for small fluctuations of the permeability, analysis based on perturbation schemes will provide reliable estimates of the effective permeability. The ensuing re-normalization methods have been applied extensively to analytically estimate the effective hydraulic conductivities of heterogeneous media (King 1988, 1996; Renard *et al.* 2000; Yeo & Zimmerman 2001; Green & Paterson 2007). Perturbation schemes for estimating effective permeability are also discussed by Drummond & Horgan (1987), Dagan (1993) and Keller (2001). A geometric averaging technique for estimating the effective permeability is presented by Jensen (1991); Abramovich & Indelman (1995) present a technique where a regularization of the Green’s function for the Laplace operator is used to calculate the effective permeability. Hristopulos & Christakos (1997) used a variational approach to determine the effective permeability for a medium where the heterogeneity exhibits a Gaussian character. The approach to upscaling adopted by Christakos (1998, 2005) is based on the concept of an epistemic cognition approach (ECA), which integrates a much wider class of information ranging from general to specific attributes. Yu & Christakos (2005) have conducted numerical experiments to generate effective permeabilities in a two-dimensional domain with anisotropic effective permeabilities.

The procedures put forward in this paper in terms of experimental techniques or data analysis are insufficiently refined to allow estimation of the full *effective symmetric* *permeability tensor* at a point corresponding to the tested cube; it is, however, instructive to calculate the *effective permeability* for the ‘cuboidal point’. Elementary measures for estimating effective permeability have been presented by a number of authors and these will be used in the calculations that follow. Several procedures for the estimation of the *effective permeability* of a porous medium are given in the literature cited previously; we employ the findings of these studies to provide estimates for the *effective permeability* (*K*_{eff}) of the cuboidal block of Indiana limestone.

Wiener (1912) bounds: the effective permeability

*K*_{eff}of a porous medium with volume*V*_{0}and permeability distribution*K*() is always bounded by the inequality 6.1 where*x**K*_{h}and*K*_{a}are, respectively, the*harmonic mean*and the*arithmetic mean*defined by 6.2Matheron (1967) suggests that the effective permeability can be estimated from the weighted average of the Wiener bounds according to 6.3 where

*α*∈[0,1] and*α*=(*D*−1)/*D*, where*D*is the dimension of the space.Journel

*et al.*(1986) proposed that the effective permeability can be derived by considering the*power average*(or an average of order*p*) with an exponent*p*in the interval between −1 and +1, depending on the spatial distribution of permeabilities, i.e. 6.4 We note that*p*=−1 corresponds to the harmonic mean, corresponds to the geometric mean and*p*=1 corresponds to the arithmetic mean. For a statistically homogeneous and isotropic medium,*p*=1−(2/*D*).The perturbation technique presented by King (1989) indicates that the effective permeability can be estimated from the result 6.5 where

*σ*^{2}is the permeability variance. Another estimate by King (1987), based on a conjecture of Landau & Lifshitz (1960) and on methods that rely on field theory, gives 6.6 where is the variance in the logarithm of permeability and*K*_{g}is the geometric mean.The estimate for the effective permeability given by Dagan (1993) takes into consideration higher order terms in the variance and has the form 6.7 where

*σ*^{2}is the variance of the permeability and*D*is the dimension of the space.

We further note that, in order to apply Matheron’s formulae, the permeability variations in the experiments should conform to a log-normal distribution. The dataset for the permeabilities given in table 1 satisfies the Kolomogorov–Smirnov test and establishes the applicability of the log-normal distribution to within a confidence level of 95%. A multiplicity of other procedures are available for estimating the effective permeability of the porous medium; for the present, however, we shall restrict attention to the estimates for the effective permeability given above and use the spatial distributions of permeability *K*(** x**) derived from the kriging exercise to evaluate the effective permeability of the cuboid of Indiana limestone. The computations were carried out at a discretized scale with average values for 64 substructured units, a typical view of which is shown in figure 13; the discretized permeabilities were obtained from the results of the kriging procedure that used the average values. The number of cubical elements used in the substructuring can be increased, but this entails an inordinate amount of data input owing to the nature of the Comsol software. The bounds proposed in equation (6.2) and the estimates defined by equations (6.3)–(6.7) can now be calculated using the values for the substructured units.

As a final calibration, we considered the substructured representation of the Indiana limestone block containing the 64 equal regions of permeability, and subjected the composite cuboidal region to one-dimensional flow in the three orthogonal directions. The three permeabilities determined from this procedure are denoted by *K*_{1},*K*_{2},*K*_{3} (these should not be interpreted as principal values of the permeability tensor). It is proposed that the *effective permeability* of the block can be estimated from the *geometric mean*
6.8

The Comsol Multiphysics finite-element code was used to perform the computations. Two opposite ends of the cuboidal region were subjected to constant potential Dirichlet boundary conditions and the remaining four sides were subjected to null Neumann conditions. Even though the externally prescribed boundary conditions correspond to a one-dimensional state of fluid flow, the internal distribution of potential is three-dimensional and will, in general, give rise to a three-dimensional flow; figure 14 clearly illustrates the development of three-dimensional flow patterns for certain choices of the hydraulic gradient. The computations provide the following estimates for ‘*apparent one-dimensional permeabilities*’ in the three orthogonal directions: , indicating a nominal measure of directional dependence. If the permeabilities in the three different directions vary by several orders of magnitude, then the interpretation of the ‘*effective permeability*’ as a geometric mean is erroneous. In the experimental results for the current research, the permeability ratio in the three orthogonal directions is 1 : 1.21 : 1.22. This is an extremely weak transverse isotropy, which suggests that the apparent mild form of transverse isotropy can be dispensed with, particularly in hydrogeological calculations, by representing the porous medium as an equivalent isotropic medium, where the effective permeability is the geometric mean. The effective estimate based on the geometric mean (6.8) gives

Table 2 presents a summary of the expressions that were used to estimate the *effective permeability*. The percentage error between the computational and the Matheron (1967) estimates is approximately 0.043 per cent and that between the computational and Journel *et al.* (1986) estimates is 0.442 per cent. Also, all estimates are within the bounds proposed by Wiener (1912). If the result of Matheron (1967) is taken as the reference, all other estimates predict the value of the *effective permeability* within an accuracy of approximately 0.45 per cent.

These estimates were also compared with the experimental results obtained previously by Selvadurai & Głowacki (2008) for Indiana limestone using axial flow tests. In its unstressed state, the permeability was estimated at (15–16)×10^{−15} m^{2}. Previous studies for Indiana limestone have also found values of permeability that range from 4×10^{−15} to 54×10^{−15} m^{2} (Churcher *et al.* 1991). The effective permeability estimates for the Indiana limestone obtained in this research point to a slightly larger value, but well within the range of variability (10^{−22} m^{2}<*K*_{eff}<10^{−14} m^{2}) that is expected of sedimentary formations consisting of limestone and dolomite (Philips 1991).

## 7. Conclusions

The determination of permeability parameters that can be used in computational models for examining subsurface fluid flow in geological media continues to pose a major challenge to many areas of environmental geosciences and geoenvironmental engineering. Theoretical and computational investigations of the annular sealing patch on a plane surface indicate that certain aspect ratios of the sealing annulus permit the use of conventional solutions applicable to a circular entry point located on the surface of a halfspace region. The experimental and theoretical results allow the mapping of the near-surface permeabilities and these can be used with a kriging procedure to determine the permeability variations within the cuboidal region. The availability of the permeability distributions of samples with larger dimensions enables a meaningful estimation of the effective permeability of the cuboidal specimen, which can provide a representative description of the point-wise isotropic effective permeability relevant to large-scale applications. It is shown that the substructured permeability distributions can be used quite conveniently to estimate several measures of effective permeability including the Wiener bounds and the estimates of Landau & Lifshitz (1960), Matheron (1967), Journel *et al.* (1986), King (1987) and Dagan (1993). For the investigations conducted, it is also found that the effective permeability of the cuboidal region of Indiana limestone can be conveniently estimated using the geometric mean. The effective permeabilities determined through these estimates confirm the validity of the Wiener (1912) bounds.

## Acknowledgements

The work described in this paper was supported in part through the Max Planck Research Prize in the Engineering Sciences awarded by the Max Planck Gesellschaft, Berlin, Germany, and in part through an NSERC Discovery Grant, both awarded to A.P.S.S. The authors are indebted to the reviewers for their helpful comments and critical reviews. The authors are also grateful to Dr V. C. Tidwell, Geohydrology Department, Sandia National Laboratories, Albuquerque, NM, USA, and to Prof. George Christakos, Department of Environmental Sciences and Engineering, University of North Carolina at Chapel Hill, NC, USA, for their advice and helpful comments.

## Footnotes

- Received September 16, 2009.
- Accepted April 14, 2010.

- © 2010 The Royal Society