3462
Ind. Eng. Chem. Res. 2002, 41, 3462-3472
Determination of the Transport Properties of Single Fractures with the Aid of Critical Path Analysis Christos D. Tsakiroglou† Foundation for Research and Technology, Institute of Chemical Engineering and High-Temperature Chemical Processes Hellas, Stadiou Street, Platani, P.O. Box 1414, GR-26500 Patras, Greece
A two-step critical path analysis (CPA) is combined with effective medium approximation for the development of phenomenological models, correlating the absolute permeability and electrical formation factor of single fractures with structural parameters of their aperture. The functional form of the CPA models is valid not only at the vicinity of percolation threshold but also far from it. The fracture aperture is regarded as a two-dimensional network of elliptical channels with the following parameters: the distributions of the major and minor channel axis lengths, the mean channel length, and the mean coordination number of the network. The models are evaluated with respect to numerical simulations in two-dimensional square lattices and experimental results from artificial glass-etched fractures and are found to be adequate for estimating the transport properties of single fractures for a broad range of parameter values of their aperture. In converse, models based on the conventional CPA, as applied to elliptical channel networks, are unable to predict satisfactorily simulated results. 1. Introduction Flow in fractures is of key importance for a variety of practical applications, such as the contaminant spreading/remediation in fractured soils, sediments and aquifers, the oil recovery from underground reservoirs, the rational exploitation of the geothermal fields situated in fractured formations, the mechanical failure in construction materials (e.g., metals and reinforced concrete), etc. Notwithstanding, it has extensively been recognized that the effective transport properties of single fractures are highly affected by their morphology; a little attention has been paid to the development of analytical phenomenological models expressing the transport properties of fractures as functions of microscopic parameters of their aperture.1 For the majority of fractured media, a single fracture rarely resembles the region bounded by two smooth parallel plates. Instead, in real fractures there may be extensive regions where the two opposite faces of the fracture surface contact each other.2 Fractured clay till is a widespread sediment type for the subsurface of North Europe.3 The analysis of scanning electron microscopy (SEM) images of resinimpregnated clay till fractures4 has revealed that fracture aperture may be regarded as a two-dimensional (2D) network of elliptical channels with parameters such as the major and minor axis length distributions, the mean channel length, and the mean coordination number. It is worthwhile to note that 2D pore networks are suitable for the description of the morphology for the majority of natural fractures encountered in soils and rocks (e.g., sandstones and limestones).1,2 Pore network models enable us to represent the structure for a variety of fractured porous media, ranging from tight impermeable to highly permeable ones, and describe transport processes with mathematically tractable models. † E-mail:
[email protected]. Fax: +3 0610 965223. Tel: +3 0610 965212.
In the past, critical path analysis (CPA)5 has successfully been used to correlate explicitly hydraulic and electrical properties of porous media with structural parameters of the pore space for very broad pore-size distributions.6-8 In addition, CPA in combination with percolation theory has been used to quantify the transport properties of heterogeneous porous media.9 Pore network simulations in 2D pore networks have been used to identify the preferential paths for hydraulic flow and electrical current in heterogeneous media and determine the geometry of the critical paths as functions of the variance of the pore-size distribution and pore connectivity.10 Moreover, CPA has been extended to calculate the permeability of fracture networks with a broad distribution of fracture aperture,11 as well as to estimate the hydraulic conductivity of a porous medium for the flow of non-Newtonian fluids.12,13 The effective medium approximation (EMA) has been found to be suitable for the determination of the electrical and hydraulic properties of “homogeneous” pore structures exhibiting narrow pore-size distributions, while CPA has been found to be complementary to EMA for the calculation of the transport properties of “heterogeneous” pore structures exhibiting high variability in pore sizes.8,12-16 Analytical relationships or numerical results, obtained with CPA and EMA, have been evaluated with respect to either experimental results or pore network simulations.7,8,10-17 Bernabe and Bruderer18 have recently evaluated the validity of various permeability models with respect to the variance of the poresize distribution and found that the flow localization occurring at high pore-size variance is not restricted to the backbone of the critical path; instead, the flow pattern is so strongly governed by the pore-size fluctuation that the use of a stochastic averaging procedure for the prediction of the transport properties of rocks from pore structure characteristics becomes inadequate. Despite the progress that has been done on the use of CPA for the correlation of the hydraulic with the electrical properties of porous media, attention has been
10.1021/ie010936w CCC: $22.00 © 2002 American Chemical Society Published on Web 06/12/2002
Ind. Eng. Chem. Res., Vol. 41, No. 14, 2002 3463
Figure 1. (a) Digitized SEM image of a 2D cross section of a resin-impregnated clay-till fracture. (b) Model of elliptical channels. (c) Major and minor channel axis length distributions of a fracture aperture as determined with image analysis.
paid mostly to simple pore shapes (e.g., cylindrical and parallel plate), whereas no analytical models have been developed to allow the independent and accurate estimation of the absolute permeability and formation factor from information concerning pore structure parameters. Such models are useful in large-scale studies of transport processes in highly heterogeneous porous media, such as, for instance, the determination of the effective transport coefficients of fractured formations through upscaling of the corresponding properties of individual fractures. In the present study, the CPA is coupled with EMA to derive analytical equations correlating the absolute permeability and electrical formation factor of single fractures with geometrical and topological parameters of their aperture. The validity of these equations is examined by comparing their predictions with computed results from pore network simulations and experimental results from artificial glass-etched single fractures. 2. Single Fractures Examination of SEM images of the cross sections of typical single fractures (Figure 1a) led us to adopt the network of elliptical channels (Figure 1b) as the most suitable model for the fracture aperture. Geometrical parameters of the fracture aperture (Figure 1a) were measured by analyzing 2D SEM images (Figure 1c) with the aid of the commercial software Sigma ScanPro.19 Topological parameters of the aperture,20 such as the mean coordination number, were estimated by identifying the skeleton of the channel network (Figure 2) and measuring the genus per node. The skeleton of the channel network of clay till single fractures was identified by combining successive 2D images of the cross section of a fracture (Figure 1a) with large-scale pictures of the fracture surface (Figure 2) and registering the channels and the interconnections between them on the picture (Figure 2). The genus of a multiply-connected surface, such as the solid/pore space interface of a
Figure 2. Skeleton of the channel network identified on the trace of a clay-till single fracture. There is a primary dense network of small channels (thin lines) intersecting with a sparse system of large root holes (thick lines). The mean network coordination number is ca. 4.0.
fracture, is a quantitative measure of its connectivity and is defined as the number of non-self-intersecting cuts that may be made upon the surface without separating it into two disconnected parts. The genus of the fracture surface is equal to the 1st Betti number, B1 of the topologically equivalent network (Figure 2) of bonds and nodes (B1 ) b - n + B0, where b is the number of bonds, n is the number of nodes, and B0 is the 0th Betti number, which coincides with the number of independent subnetworks).4,20 The aperture of a single fracture is modeled as a disordered 2D network of elliptical channels (Figure 1b) of equal length, lp. The major, l1, and minor, l2, axis lengths are assigned to network channels completely at random and according to two different distribution functions, f1(l1) and f2(l2), respectively. The channel connectivity can be quantified by the mean coordination number, z, and lattice type. For spatially uncorrelated pore networks,21 these two topological parameters suffice to specify the bond percolation threshold, pc, which
3464
Ind. Eng. Chem. Res., Vol. 41, No. 14, 2002 Table 1. Structural Properties of Artificial Fractures
Figure 3. Patterns of channel network models.
property
pattern M-1
pattern M-2
topology 〈l1〉 (µm) σl1 (µm) 〈l2〉 (µm) σl2 (µm) lp (µm)
square lattice 470 168 126 20 1365
square lattice 774 274 125 20 1840
(Table 1) are comparable to the corresponding ones of real fractures. However, the variability of the minor and major channel axis lengths of artificial fractures is smaller than that of real ones (Table 1 and Figure 1c). The volumetric flow rate, Q, of a liquid (oil, water, or brine) injected in the artificial fracture is adjusted by using a syringe pump (HARVARD PHD2000), while the pressure drop, ∆P, is measured along a length L of the central region of the network by using a differential pressure transducer (VALIDYNE DPD V/P55D). The absolute permeability, k, is estimated by fitting the experimental measurements Q - ∆P to Darcy’s law, given by23
Q)
Ak ∆P µL
(1)
where µ is the liquid viscosity, A is the total (apparent) cross-sectional area of the fracture, defined by
A ) Wlp
(2)
W is the total fracture width (Figure 3), and lp is the length of periodicity (Figures 1b and 4), namely, the distance between the centers of two adjacent nodes. The formation factor of each artificial fracture was determined by saturating the channel network with several aqueous solutions of NaCl in different concentrations and measuring, with the use of two miniaturized electrodes, the electrical conductivity, c0, along a thin strip of length L and cross-sectional area AF (which, in practice, coincides with the mean pore cross-sectional area, namely, AF ) (π/4)〈l1〉〈l2〉), as a function of the solution conductivity, cw. With reference to the apparent pore cross-sectional area, Ak ) lp2, the formation factor is defined by
F)
( )( ) Ak c0 AF c w
(3)
and calculated from the slope of the c0-cw relationship. Figure 4. Detailed morphology of artificial fractures fabricated with the pattern (a) M-1 (lp ) 1365 µm) and (b) M-2 (lp ) 1840 µm).
is defined as the smallest fraction of large bonds (channels) which ensure the creation of a noninterrupting pathway (network-spanning cluster) between the inlet and outlet of the network.15,22 Two artificial large 2D single fractures (dimensions ∼ 15 cm × 11 cm) of controlled morphology were fabricated by etching mirror image patterns (Figure 3) of channel networks on two glass plates with photolithographic techniques and sintering the glass plates in a programmable furnace23 (Figure 4). The channels have a lenticular cross-sectional shape, which is very similar to the elliptical one.24 The channel dimensions and coordination number (z ) 4) of artificial fractures
3. CPA The fluid flow or electrical conduction through a pore system with broad pore-size distribution is governed by a critical (hydraulic or electrical) pore conductance, gc, which is defined by the condition that all pores with conductance greater than gc create a network-spanning cluster transferring most of the flow or current.6-8 This critical conductance gc depends on the structural characteristics (e.g., pore shape, pore-size distribution, connectivity, dimensionality, etc.) of the pore space, and CPA is used to calculate the corresponding critical pore dimensions. According to CPA,6 all conductances with values greater than gc are set equal to gc, while all conductances smaller than gc are set equal to zero. The transport properties of the critical path may be dominated by universal scaling laws prevailing in the vicinity
Ind. Eng. Chem. Res., Vol. 41, No. 14, 2002 3465
of the percolation threshold.12,15,22,25 In mathematical terms, the fluid flow or electrical conduction problem reduces to the calculation of the critical pore dimensions that maximize the total conductance of the critical path.6 In general, the dependence of the local hydraulic conductance on pore size may be different from that of the local electrical conductance, so that different critical pore dimensions may be necessary to be determined for the two transport problems.6-10 With reference to the aperture model of real fractures, the hydraulic conductance for laminar flow in elliptical channels, gh, is given by
l13l23
gh ) ch 2 l1 + l22
(lh1c)3l2c3
(4)
where
ch ) π/64µlp
(5)
The electrical conductance of elliptical channels, ge, is given by
ge ) cel1l2
(7)
∫l∞f1(l1) dl1 ) pc
(8a)
∫l∞f2(l2) dl2 ) pc
(8b)
2c
where pc is the network percolation threshold, which, for bond percolation problems in 2D systems is a function of network coordination number, and can be approximated by the relation15,25
pc ) 2/z
(9)
while details of the network topology (lattice type) are reflected in the corresponding accessibility function.21,25 In other words, the lengths l1c and l2c correspond to the maximum values of channel axes that ensure the creation of two different and uncorrelated critical paths of minimum connectivity (percolating clusters) throughout the pore network.5,8 Given that the variance of the minor axis length, σl2/〈l2〉, is usually smaller than that of the major axis length, σl1/〈l1〉 (Figure 1c and Table 1),
(10)
then a conductance equal to ghc is assigned to channels with l1 g lh1c and l2 g l2c, whereas a conductance equal to zero is assigned to channels with l1 < lh1c and l2 < l2c. At the vicinity of the percolation threshold, the overall conductance of the hydraulic critical path, σh, follows a universal scaling law22,25 and may be written as
σh ) ghc
The straightforward application of the classical CPA to a network of elliptical pores includes the determination of two pairs of critical lengths for flow (lh1c and lh2c) and conduction (le1c and le2c) so that the overall hydraulic and electrical conductances of the corresponding critical paths are maximized.6,7 This conventional approach is complicated by the dependence of gh and ge on two variables, l1 and l2, and is presented in detail at the end of this section. Alternatively, the problem can be simplified significantly if CPA is applied separately with respect to each length. According to percolation theory, two critical (percolation) channel axis lengths, l1c and l2c, can be defined so that channels with l1 > l1c and l2 > 0 or l2 > l2c and l1 > 0 create a continuous networkspanning cluster for a first time, namely
1c
ghc ) ch h 2 (l1c) + l2c2
(6)
where
ce ) πcw/4lp
it is assumed that the critical minor axis length is equal to l2c for both problems of liquid flow and electrical conduction. In light of the foregoing simplification, the CPA problem is reduced to the determination of the critical major channel axis lengths, lh1c e l1c and le1c e l1c, that maximize the hydraulic and electrical conductance of connected pathways of channels having l1 g lh1c and l1 g le1c, respectively. By considering that the fluid flow is controlled by a critical hydraulic conductance ghc given by
[
]
t
p1(lh1c) - pc 1 - pc
(11)
where t is the universal scaling exponent of conductivity (t = 1.3 for 2D systems22), p1(l1) denotes the probability that a channel has a major axis length greater than l1 and is given by
p1(l1) )
∫l∞f1(l1′) dl1′ 1
(12)
and the term 1 - pc is a normalization factor ensuring that the overall conductance tends to 0 and ghc at the limits p1 f pc and p1 f 1, respectively. Likewise, assuming that the electrical conduction in the network is controlled by a critical electrical conductance, gec, given by
gec ) cele1cl2c
(13)
the overall electrical conductance of the channel network can be approximated by the relation
σe ) gec
[
]
p1(le1c) - pc 1 - pc
t
(14)
Note that σh and σe are microscopic quantities and should be regarded as the effective values of the hydraulic and electrical conductance of the critical path per each channel. Far from the vicinity of percolation threshold (p1 . pc), the microscopic effective conductance of a connected pathway of channels with length greater than a critical value can be deduced by using the EMA. According to EMA, the individual conductances of pores are replaced by a mean effective value supposing that the average fluctuation caused on the (potential or pressure) field by this replacement is equal to zero.15,16,25 Having set in CPA that the network consists of a fraction of channels p1 with conductance ghc or gec (hydraulic or electrical critical path) and a fraction of channels
3466
Ind. Eng. Chem. Res., Vol. 41, No. 14, 2002
1 - p1 with zero conductance, the EMA finally yields the microscopic mean effective conductance25
zp1(lh1c) - 2 σh ) ghc z-2
(15)
for hydraulic flow and
zp1(le1c) - 2 gec σe ) z-2
[ [
σh ) ghc
σe ) gec
] ]
- pc 1 - pc
p1(le1c)
- pc 1 - pc
dlh1c
)0
dσe dle1c
)0
(18)
∫
dP dl
cwlp2 dP im ) F dl
(27)
respectively. It is well-known that the flow (or current) transferred by the critical path is a substantial or small portion of the total flow (or current), depending on the specific geometry and topology of the investigated pore system.10,17,18 To be able to derive analytical expressions of the absolute permeability and formation factor as functions of the critical channel axis lengths, we adopt the approximate assumption that the ratio of flow rates QCPA/Qm and currents iCPA/im are roughly equal to the corresponding ratio of flow areas, namely,
(28)
(29)
(19)
iCPA AeCPA p1(le1c) le1c = = im Am 〈l1〉
(20)
respectively. The foregoing assumption is equivalent to the condition that the mean pore flow velocity and mean pore current flux in a disordered network are approximately equal to the corresponding ones prevailing in the hydraulic and electrical critical paths. By combining eqs 22-29 and taking into account eqs 11 and 14, finally we obtain the approximate relationships
(21)
Equations 20 and 21 are easily solved numerically with a Newton-Raphson method so that the critical lengths lh1c and le1c are determined. Under a pressure gradient dP/dl, the microscopic volumetric flow rate (namely, flow per channel) across the hydraulic critical path, QCPA, is given by
QCPA ) -σh′
(26)
and
and
k=
1c
klp2 dP µ dl
QCPA AhCPA p1(lh1c) lh1c = = Qm Am 〈l1〉
h 1c
e e e f (l ) dl1 - pc) - tl1c f1(l1c) ) 0 l 1 1
(25)
In addition, the microscopic volumetric rate, Qm, and current, im, for the entire network are defined by
Qm ) -
∫l∞f1(l1) dl1 - pc) - tlh1c f1(lh1c) ) 0
∞
(24)
σe′ ) σelp
(17)
and
(
dV dl
where
which in combination with eqs 10-14 result in
(lh1c)2 + 3l2c 2 ( (lh1c)2 + l2c2
(23)
Also, under a voltage gradient dV/dl, the microscopic current (namely, current per channel) across the electrical critical path, iCPA, is given by
(16)
which are completely analogous to eqs 11 and 14, by using the exponent value t ) 1.0. Consequently, it becomes clear that, for 2D lattices, the functional form of CPA equations is valid not only at the vicinity of percolation threshold but also far from it, with the difference that the exponent specifying the dependence of the microscopic effective conductance on the fraction of channels of the critical path is not constant but varies from its universal value, t ) 1.3 (p1 f pc) to t ) 1 (p1 f 1). In other words, the exponent t is a decreasing function of the fraction of channels of the critical path, p1 (dt/dp < 0). As lh1c f l1c and le1c f l1c, then p1(lh1c) f pc and p1(le1c) f pc and hence σh f 0 and σe f 0. As lh1c f 0 and le1c f 0, then ghc f 0 and gec f 0 and hence σh f 0 and σe f 0. Therefore, for some intermediate values of the critical lengths, 0 < lh1c < l1c and 0 < le1c < l1c, the conductances σh and σe are maximized. The maximization of σh and σe with respect to lh1c and le1c is based on the conditions
dσh
σh′ ) σhlp
iCPA ) -σe′
for electrical conduction, where z is the network coordination number. In conjunction with eq 9, eqs 15 and 16 are equivalently written as
p1(lh1c)
where
(22)
(
〈l1〉
p1(lh1c) lh1c F=
(
)( )(
)[ )( )( )[
]
(lh1c)3l2c3 p1(lh1c) - pc π 1 - pc 64lp2 (lh1c)2 + l2c2
p1(le1c) le1c 4lp2 1 - pc 1 e π l l p (le ) - p 〈l1〉 1c 2c 1 1c c
]
t
(30)
t
(31)
for the network absolute permeability and formation factor, respectively. In general, the absolute permeability of a network of pores of arbitrary cross-sectional shape is scaled as
k ∝ rh,eff4/lp2
(32)
Ind. Eng. Chem. Res., Vol. 41, No. 14, 2002 3467
where rh,eff is an effective value of the equivalent hydraulic radius, rh, defined by
rh ) 4(Ap/Pp)
(33)
and Ap and Pp are the area and perimeter of a pore. Obviously, the radius rh,eff can be determined by using either EMA for narrow distributions of pore dimensions or CPA for wide distributions of pore dimensions.8,10,16 For elliptical channels with l1 . l2, the hydraulic radius rh,eff is a strong function of the minor axis length and a weak function of the minor to major axis length, and such a behavior is completely consistent with eq 30, which implies that
k∝
〈l1〉l2c
The application of the conventional CPA approach to elliptical channel networks results in the following relationships for the critical hydraulic and electrical conductance
(34)
Moreover, the formation factor of a pore network is proportional to
F ∝ lp2/re,eff2
(35)
where re,eff is an effective value of the equivalent electrical radius, re, defined by
re ) (Ap/π)
1/2
(36)
Evidently, for elliptical channels, the radius re is a strong function of the minor and major axis lengths, while eq 31 yields
F ∝ lp2/〈l1〉l2c
(37)
This means that eq 31 is expected to fail at increasing high σl1/〈l1〉 values, because it is unable to account for the effects of l1c and le1c on F values. Alternatively, the analysis of the electrical critical path can be repeated by assuming a critical major channel axis length equal to l1c, determining the critical minor axis length, le2c, from the solution of equation
(
∫l f2(l2) dl2 - pc) - tle2cf2(le2c) ) 0 ∞
and finally obtaining the relationship
F=
(38)
e 2c
( )( )( )[
p(le2c)le2c 4lp2 1 - pc 1 e π 〈l2〉 l2cl1c p2(le2c) - pc
]
p2(l2) )
∫l f2(l2′) dl2′ 2
gec ) cele1c le2c
(43)
[ [
] ]
p(lh1c,lh2c) - pc 1 - pc
t
σh ) ghc
p(le1c,le2c) - pc 1 - pc
t
σe ) gec
(44)
(45)
In the foregoing relations, p(l1,l2) denotes the probability that a channel has major and minor axis lengths greater than l1 and l2, respectively, and is defined by
p(l1,l2) )
∫l∞∫l∞f(l1′,l2′) dl1′ dl2′ 2
1
(46)
where f(l1,l2) is a joint distribution function that reduces to
f(l1,l2) ) f1(l1) f2(l2)
(47)
by ignoring any correlation of l1 to l2. The maximization of the overall conductance of the hydraulic and electrical critical path is based on the conditions
∂σh ∂lh1c
)0
∂σh ∂lh2c
)0
(48)
)0
(49)
and
∂σe ∂le1c
)0
∂σe ∂le2c
respectively. Equation 48 in conjunction with eqs 42, 44, 46, and 47 results in
t
(39)
where p2(l2) denotes the probability that a channel has a minor axis length greater than l2 and is given by ∞
(42)
and overall conductance of the hydraulic and electrical critical path
3
[1 + (l2c/lh1c)2]lp2
(lh1c)3(lh2c)3 ghc ) ch h 2 (l1c) + (lh2c)2
(40)
(lh1c)2 + 3(lh1c)2 [p1(lh1c) p2(lh2c) - pc] h 2 h 2 (l1c) + (l2c) tlh1cf1(lh1c) p2(lh2c) ) 0 (50a) (lh2c)2 + 3(lh2c)2 [p1(lh1c) p2(lh2c) - pc] h 2 h 2 (l1c) + (l2c) tlh2cf2(lh2c) p1(lh1c) ) 0 (50b)
Equation 39 implies that
(41)
Accordingly, eq 49 in conjunction with eqs 43, 45, 46, and 47 results in
and hence the estimated F values are expected to be sensitive to the variance, σl1/〈l1〉 (namely, to changes of l1c) at increasing values of this parameter. Therefore, eqs 31 and 39 seem to be adequate for narrow and wide distributions of the major channel axis length, respectively.
p1(le1c) p2(le2c) - pc - tle1cf1(le1c) p2(le2c) ) 0 (51a)
F ∝ lp2/〈l2〉l1c
p1(le1c) p2(le2c) - pc - tle2cf2(le2c) p1(le1c) ) 0 (51b) Numerical solution of the nonlinear systems of eqs 50
3468
Ind. Eng. Chem. Res., Vol. 41, No. 14, 2002
and 51 provides the critical lengths (lh1c, lh2c) and (le1c, le2c), which, in turn, can be used for the calculation of the absolute permeability and formation factor through the relationships
k=
(
〈l1〉〈l2〉
p1(lh1c) p2(lh2c) lh1c lh2c
)( )( [
)
(lh1c)3(lh2c)3 π × 64lp2 (lh1c)2 + (lh2c)2
]
p1(lh1c) p2(lh2c) - pc 1 - pc
and
F=
(
)( )( ) [
p1(le1c) p2(le2c) le1c le2c 4lp2 1 × π le le 〈l1〉〈l2〉 1c 2c 1 - pc
p1(le1c) p2(le2c) - pc
t
(52)
]
t
(53)
respectively. 4. Pore Network Simulation The predictability of the foregoing models can easily be evaluated with respect to k and F values obtained with simulations on pore networks. A major axis length, l1, and a minor axis length, l2, are randomly chosen according to the predefined frequency distribution functions, f1(l1) and f2(l2), and are randomly assigned to each channel of a regular square lattice. When a macroscopic voltage is imposed along the network, the potentials at nodes are calculated by using Kirchoff’s conservation law, formulating that the algebraic sum of the currents flowing into each node is equal to zero. From the application of this rule to each node of the network, a system of coupled linear equations for the potentials at nodes is obtained.7,12,13,16,25 The system is solved with a properly modified Gauss-Seidel numerical method so that the total current flow at the network inlet/outlet and therefrom the network electrical conductance and formation factor are calculated. The problem of hydraulic flow is solved by using exactly the same method, the analogies pressure-potential, flow rate-current, hydraulic conductance-electrical conductance, and calculating the total flow rate at the network inlet/outlet and therefrom the absolute permeability. Because of the statistical nature of theoretical simulations, a sufficient number of realizations and a large network are commonly required in order to minimize the uncertainties commonly embedded in calculated k and F values. 5. Results and Discussion A sensitivity analysis in 2D square lattices with lognormal distributions f1(l1) and f2(l2) revealed that the critical lengths l1c, lh1c, and le1c and critical probabilities p(lh1c) and p(le1c) are decreasing functions of the variance σl1/〈l1〉 (Figures 5 and 6). According to eqs 10 and 13, the hydraulic effective conductance, σh, is affected more strongly than the electrical effective conductance, σe, by the critical major channel axis length, and unavoidably it takes on a local maximum at a length value of lh1c g le1c (Figure 5a,b). The nonlinear dependence of the critical major axis lengths on the variance σl1/〈l1〉, within a broad value range of this parameter (Figure 5), may
Figure 5. Sensitivity of the critical major axis lengths to the variance σl1/〈l1〉: (a) σl2/〈l2〉 ) 1/3; (b) σl2/〈l2〉 ) 1.0.
be reflected in an analogous behavior of k and F. Therefore, it would be interesting to investigate the capability of models to predict the transport properties of single fractures within a broad range of values of the variances σl2/〈l2〉 and σl1/〈l1〉. The electrical critical path has a higher connectivity than that of the percolating cluster, over the entire region of σl1/〈l1〉 values, and tends asymptotically to the percolation cluster [p1(le1c) f pc] at very high σl1/〈l1〉 values (Figure 6). Although the connectivity of the hydraulic critical path is also higher than that of the percolating cluster, over the widest region of σl1/〈l1〉 values, this critical path has the tendency to “enter” into the vicinity of the percolation threshold [p1(lh1c) f pc] at high but finite σl1/〈l1〉 values (Figure 6). Consequently, the best estimates of F and k are expected for t ) 1.0 or 1.3, depending on whether the variance σl1/〈l1〉 is low or high, respectively. Simulations were performed in square channel networks for the determination of accurate k and F values that might be used for the evaluation of the CPA models. The size of the network is associated with the uncertainty embedded into k and F values but does not have any pronounced effect on their average values (Table 2). In the simulations, networks of dimensions 60 × 30 with log-normal distribution functions f1(l1) and f2(l2) and various σl1/〈l1〉 and σl2/〈l2〉 values were used. In each case, the mean value and standard deviation of k and 1/F were computed by averaging over the results of 15
Ind. Eng. Chem. Res., Vol. 41, No. 14, 2002 3469
Figure 6. Sensitivity of the critical probabilities to the variance σl1/〈l1〉: (a) σl2/〈l2〉 ) 1/3; (b) σl2/〈l2〉 ) 1.0. Table 2. Effect of Network Size on Simulated k and F Values 〈l1〉 ) 1000 µm, σl1/〈l1〉 ) 0.5, 〈l2〉 ) 300 µm, and σl2/〈l2〉 ) 0.5 network size
k (Da)
F
30 × 15 60 × 30 100 × 50
324.3 ( 26.0 326.6 ( 10.2 324.1 ( 4.85
11.97 ( 0.46 11.86 ( 0.17 11.94 ( 0.11
realizations (Figure 7a,b). The inverse formation factor, 1/F, shown in Figure 7b is a quantitative measure of the electrical conductance of the pore network and can be compared directly to permeability, k (Figure 7a), which is a measure of the hydraulic conductance of the pore network. Simulated knet and Fnet values were compared with corresponding model estimates, kmod and Fmod, in terms of the relative errors defined by
Ek )
|
knet - kmod knet
|
(54)
EF )
|
Fnet - Fmod Fnet
|
(55)
and
respectively, and the results are shown in Table 3. Consistent with the earlier described behavior of critical paths, the exponent value t ) 1.0 seems to be the best
Figure 7. Comparison of predictions of models resulting from the two-step CPA (eqs 30, 31, and 39) with pore network simulations for the (a) absolute permeability and (b) formation factor.
choice for estimating the absolute permeability of networks with low and medium σl1/〈l1〉 and σl2/〈l2〉 values, while the exponent value t ) 1.3 is adequate for estimating the absolute permeability of networks with high σl1/〈l1〉 and σl2/〈l2〉 values (Table 3 and Figure 7a). In addition, the model of eq 31 is satisfactory for quite low values of the variance σl1/〈l1〉 and t ) 1.0 but fails to produce reliable estimates of the formation factor for pore networks with high σl1/〈l1〉 values (Table 3 and Figure 7b). In contrast, the predictions of eq 39 agree satisfactorily with pore network simulations at high σl1/ 〈l1〉 values and are improved only slightly if the exponent value t ) 1.3 is used (Table 3 and Figure 7b). Finally, estimates of the models of eq 30 and eqs 31 and 39 are compared with the experimental F and k values of two artificial fractures (Table 4) with Gaussian distributions of major and minor channel axis lengths (Table 1). Because of the relatively small variance of the major and minor channel axis lengths for both fractures (Table 1), both models of eqs 31 and 39 produce satisfactory estimates of F independently on t value, whereas the best estimates of k, eq 30, are taken for t ) 1.0 (Table 4). Models based on the conventional CPA, eqs 50-53, are unable to predict satisfactorily the transport properties of 2D networks of elliptical channels (Figure 8a,b). In practice, these models subestimate the values of the critical channel axis lengths or overestimate the critical
3470
Ind. Eng. Chem. Res., Vol. 41, No. 14, 2002
Table 3. Deviation of Model Predictions (Eqs 30, 31, and 39) from Pore Network Simulations (lp ) 1500 µm) σl2/〈l2〉 ) 1/3 Ek (%) (eq 30)
EF (%) (eq 31)
σl2/〈l2〉 ) 1.0 EF (%) (eq 39)
Ek (%) (eq 30)
EF (%) (eq 31)
EF (%) (eq 39)
σl1/〈l1〉
t)1
t ) 1.3
t)1
t ) 1.3
t)1
t ) 1.3
t)1
t ) 1.3
t)1
t ) 1.3
t)1
t ) 1.3
0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.7 2.0 2.5
10.5 14.6 15.5 17.9 14.6 12.5 10.7 8.1 2.1 0.7
11.1 18.5 22.2 27.2 26.8 27.3 27.9 28.6 26.4 27.5
2.9 3.0 0.8 3.3 9.4 14.7 19.4 26.9 33.8 42.1
4.0 4.6 1.1 0.4 6.3 11.1 15.6 22.9 29.8 38.1
3.6 3.3 1.6 3.6 2.4 2.5 2.9 2.7 1.6 2.4
4.9 4.6 2.8 5.1 3.8 3.8 4.4 4.0 2.9 3.7
3.6 3.6 4.8 1.3 7.9 13.5 25.1 31.4 35.5 62.2
4.9 7.0 10.7 7.8 4.8 2.7 4.3 5.5 5.4 20.3
1.7 1.7 0.6 4.5 10.4 15.7 21.6 27.4 32.0 42.5
2.5 3.4 1.2 1.7 7.1 12.1 18.0 23.4 28.0 38.6
9.8 10.0 9.9 9.5 10.4 10.5 11.5 9.9 7.7 10.3
6.6 6.7 6.6 6.3 7.2 7.2 8.4 6.6 4.5 7.1
Table 4. Comparative Results of Model Predictions and Experimental Results pattern kexp, Da Fexp kmod, Da (Ek) Fmod (EF) (eq 31) Fmod (EF) (eq 39)
t ) 1.0 t ) 1.3 t ) 1.0 t ) 1.3 t ) 1.0 t ) 1.3
M-1
M-2
20.26 44.42 17.15 (15.3%) 15.96 (21.2%) 47.26 (6.4%) 49.15 (10.6%) 42.13 (5.2%) 45.96 (3.5%)
19.34 50.16 17.13 (11.4%) 16.24 (16.0%) 52.46 (4.6%) 54.13 (7.9%) 46.9 (6.5%) 47.4 (5.5%)
probabilities because they are based on the restrictive assumption that the hydraulic (or electrical) critical path consists of a sufficient fraction of channels with l1 g lh1c and l2 g lh2c (or l1 g le1c and l2 g le2c) or equivalently p1(lh1c) p(lh2c) g pc [or p1(le1c) p(le2c) g pc]. The phenomenological models derived here can bridge the gap between the microscopic characterization of fractured media, which is carried out with the use of laboratory-scale analyses at the scale of a single fracture, and the geological characterization of fractured formations, which is done with field work and the use of geophysical techniques at the scale of fracture networks. Specifically, the present models can be used as fast and cost-effective tools for the introduction of information concerning the conductivity of single fractures into an upscaling scheme properly designed to estimate the large-scale transport properties (e.g., absolute and relative permeabilities, formation factor, resistivity index, etc.) of fractured rocks and soils. 6. Conclusions A 2D disordered network of elliptical channels is considered as a fracture aperture model to develop approximate analytical models for the absolute permeability and formation factor of single fractures. Ideas from percolation theory and EMA are used in the development of a two-step CPA. The new models are evaluated with respect to pore network simulations and experimental results from glass-etched artificial fractures. The main conclusions are outlined below: (i) The models are capable of producing reliable estimates for k and F within a broad spectrum of geometrical parameters of the fracture aperture and specifically within a broad range for the variance of the major and minor channel axis length distributions. (ii) The electrical and hydraulic critical paths are usually far from the neighborhood of the percolation threshold, and their effective conductivity is governed by EMA (t ) 1.0) rather than universal scaling. At increasing values of the variance of channel axis lengths, the critical paths gradually get near the percolation
Figure 8. Comparison of predictions of models resulting from the conventional CPA (eqs 52 and 53) with pore network simulations for the (a) absolute permeability and (b) formation factor.
threshold, and their effective conductance is dominated by the universal scaling exponent of conductivity (t ) 1.3). (iii) Despite their approximate character, the new models agree satisfactorily with theoretical and experimental results. (iv) Models derived from the application of the conventional CPA to networks of elliptical channels are unable to reproduce simulated results. (v) The models are useful tools in the determination of the transport properties of single fractures from microscopic characteristics of their aperture and hence
Ind. Eng. Chem. Res., Vol. 41, No. 14, 2002 3471
in upscaling of the transport properties of heterogeneous fracture networks. Acknowledgment This work was performed under Energy Environment and Sustainable Development (EESD) Contract No. EVK1-CT1999-00013 (project acronym: TRACe-Fracture) supported by the European Commission. Thanks are due to my colleagues, Dr. M. Theodoropoulou and Dr. V. Karoutsos, for their assistance in performing hydraulic and electrical tests on artificial fractures. Mr. K. E. Klint from the Geological Survey of Denmark and Greenland (GEUS) is also acknowledged for the provision of SEM images and identification of the channel network skeleton on clay-till single fracture. Nomenclature A ) apparent cross-sectional area of fracture, m2 h e ACPA (ACPA ) ) cross-sectional area per channel of the hydraulic (electrical) critical path, m2 AF ) mean channel cross-sectional area, m2 Ak ) apparent channel cross-sectional area, m2 Am ) mean cross-sectional area per channel of network with l2 ) l2c, m2 Ap ) individual pore cross-sectional area, m2 c0 ) electrical conductivity of fracture completely saturated with electrolyte, S/m cw ) electrolyte conductivity, S/m f(l1,l2) ) joint distribution function of major and minor channel axis lengths, m-2 f1(l1) ) major channel axis length distribution function, m-1 f2(l2) ) minor channel axis length distribution function, m-1 F ) formation factor gc ) critical (hydraulic or electrical) conductance, S (Ω-1) ge ) electrical conductance, S (Ω-1) gec ) critical electrical conductance, S (Ω-1) gh ) hydraulic conductance, m3 s-1 Pa-1 ghc ) critical hydraulic conductance, m3 s-1 Pa-1 iCPA ) mean current per channel of the critical path, A im ) mean current per channel of the network, A k ) absolute permeability, Da (1 Da ) 0.987 × 10-12 m2) lp ) channel length, m l1 (l2) ) major (minor) channel axis length, m l1c (l2c) ) percolation length of major (minor) channel axis, m le1c (le2c) ) critical length of major (minor) channel axis for the electrical path, m lh1c (lh2c) ) critical length of major (minor) channel axis for the hydraulic path, m 〈l1〉 (〈l2〉) ) mean value of the major (minor) channel axis length, m L ) macroscopic fracture length, m p ) percolation probability for both channel axis lengths pc ) percolation threshold p1 ) percolation probability for the major axis length p2 ) percolation probability for the minor axis length Pp ) pore perimeter, m Q ) volumetric flow rate, m3 s-1 QCPA ) mean flow rate per channel in the hydraulic critical path, m3 s-1 Qm ) mean flow rate per channel in the network, m3 s-1 rh (re) ) equivalent hydraulic (electrical) pore radius, m t ) universal scaling exponent of conductivity W ) fracture width, m z ) network coordination number (number of bonds per node)
Greek Letters ∆P ) pressure drop, Pa µ ) liquid viscosity, kg m-1 s-1 σe ) effective conductance of the electric critical path, S σh ) effective conductance of the hydraulic critical path, m3 s-1 Pa-1 σl1 (σl2) ) standard deviation of the major (minor) channel axis length, m Subscripts eff ) effective value exp ) experimentally measured value mod ) value calculated with CPA model net ) value calculated with network simulation
Literature Cited (1) Adler, P. M.; Thovert, J.-F. Fracture and Fracture Networks; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1999. (2) Zimmerman, R. W.; Bodvarsson, G. S. Hydraulic Conductivity of Rock Fractures. Transp. Porous Media 1996, 23, 1. (3) Nilsson, B.; Sidle, R. C.; Klint, K. E.; Boggild, C. E.; Broholm, K. Mass Transport and Scale-Dependent Hydraulic Tests in a Heterogeneous Glacial Till-Sandy Aquifer System. J. Hydrol. 2001, 243, 162. (4) Klint, K. E.; Tsakiroglou, C. D. A New Method of Fracture Aperture Characterization. Proceedings of the 5th International Conference on Protection and Restoration of the Environment, Thassos Island, Greece, July 2-6, 2000; Vol. 1, p 127. (5) Ambegaokar, V.; Halperin, N. I.; Langer, J. S. Hopping Conductivity in Disordered Systems. Phys. Rev. B: Solid State 1971, 4, 2612. (6) Katz, A. J.; Thompson; A. H. Quantitative Prediction of Permeability in Porous Rock. Phys. Rev. B: Condens. Matter Mater. Phys. 1986, 34, 8179. (7) Katz, A. J.; Thompson, A. H. Prediction of Rock Electrical Conductivity from Mercury Injection Measurements. J. Geophys. Rev. 1987, B92, 599. (8) Friedman, S. P.; Seaton, N. A. Critical Path Analysis of the Relationship between Permeability and Electrical Conductivity of Three-Dimensional Pore Networks. Water Resour. Res. 1998, 34 (7), 1703. (9) Hunt, A. G. Applications of Percolation Theory to Porous Media with Distributed Local Conductances. Adv. Water Resour. 2001, 24, 279. (10) David, C. Geometry of Flow Paths for Fluid Transport in Rocks. J. Geophys. Rev. 1993, 98, 12267. (11) Charlaix, E.; Guyon, E.; Roux, S. Permeability of a Random Array of Fractures of Widely Varying Apertures. Transp. Porous Media 1987, 2, 31. (12) Sahimi, M. Nonlinear transport processes in disordered media. AIChE J. 1993, 39, 369. (13) Shah, C. B.; Yortsos, Y. C. Aspects of flow of power-law fluids in porous media. AIChE J. 1995, 41, 1099. (14) Friedman, S. P.; Seaton, N. A. On the Transport Properties of Anisotropic Networks of Capillaries. Water Resour. Res. 1996, 32, 339. (15) Sahimi M. Flow and Transport in Porous Media and Fractured Rock; VCH: Weinheim, Germany, 1999. (16) David, C.; Gueguen, Y.; Pampoukis, G. Effective Medium Theory and Network Theory Applied to the Transport Properties of Rock. J. Geophys. Rev. 1990, 95, 6993. (17) Bernabe, Y. The Transport Properties of Networks of Cracks and Pores. J. Geophys. Res. 1995, 100, 4231. (18) Bernabe, Y.; Bruderer, C. Effect of the Variance of Pore Size Distribution on the Transport Properties of Hetreogeneous Media. J. Geophys. Res. 1998, 103, 513. (19) Klint, K. E.; Rosenboom, A.; Tsakiroglou, C. D. A New Approach to Analysis of Mechanical Fracture Aperture. In Fractured Rock 2001; Kueper, B. H., et al., Eds.; Proceedings of the International Conference “Fractured Rock 2001”, Toronto, Canada, 2001; Queen’s University: Kingson, Ontario, Canada, 2001. (20) Tsakiroglou, C. D.; Payatakes, A. C. Characterization of the Pore Structure of Reservoir Rocks with the aid of Serial Sectioning Analysis, Mercury Porosimetry and Network Simulation. Adv. Water Resour. 2000, 23, 773.
3472
Ind. Eng. Chem. Res., Vol. 41, No. 14, 2002
(21) Yanuka, M. Percolation Processes and Porous Media. II. Computer Calculations of Percolation Probabilities and Cluster Formation. J. Colloid Interface Sci. 1989, 127, 35. (22) Stauffer, D.; Aharony, A. Introduction to Percolation Theory, 2nd ed.; Taylor & Francis: London, 1992. (23) Theodoropoulou, M.; Karoutsos, V.; Tsakiroglou, C. D. Investigation of the Contamination of Fractured Formations by Non-Newtonian Oil Pollutants. Environ. Forensics 2001, 2, 321334. (24) Tsakiroglou, C. D.; Payatakes, A. C. Mercury Intrusion and Retraction in Model Porous Media. Adv. Colloid Interface Sci. 1998, 75, 215.
(25) Tsakiroglou, C. D.; Fleury, M. Pore network analysis of resistivity index for water-wet porous media. Transp. Porous Media 1999, 35, 89. (26) Kozicki, W.; Chou, C. H.; Tiu, C. Non-Newtonian flow in ducts of arbitrary cross-sectional shape. Chem. Eng. Sci. 1966, 21, 665.
Received for review November 19, 2001 Revised manuscript received April 8, 2002 Accepted April 15, 2002 IE010936W
