A THEORETICAL TREATMENT OF MICROSCOPIC PHENOMENA IN POROUS ROCKS by ANTHONY LEE ENDRES A.Sc. Science, Grand Rapids Junior College, 1975 B.Sc.(High Honors) Mathematics, Michigan Technological University, 1977 M.Sc. Mathematics, Texas A&M University, 1979 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF GEOPHYSICS AND ASTRONOMY THE INSTITUTE OF APPLIED MATHEMATICS AND STATISTICS We accept this thesis as conforming to the required standards THE UNIVERSITY OF BRITISH COLUMBIA December 1991 Â© Anthony Lee Endres In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. (Signature) Department of Geophysics and, Astronomy The University of British Columbia Vancouver, Canada Date December 5, 1991 DE.6 (2/88) Abstract This thesis makes five significant contributions to the theoretical analysis of the connection between microscopic phenomena and the macroscopic physical properties of porous rocks. A random sphere packing model permitting contact generation during hydrostatic compression is derived. It is demonstrated that the closure of near-contact gaps with an extremely small mean width significantly alters the elastic properties of granular media from those predicted by previous models in the pressure range used in laboratory measurements. Generalized forms of the inclusion-based formulations are obtained; major classes of these formulations are defined according to the manner in which interactions between inclusions in the heterogeneous system are simulated. Each class possesses an associated microstructure that determines the topological relationship between the various components. Inclusion-based formulations are obtained that describe the effects of pore-scale fluid distribution on the dielectric and elastic responses of a partially saturated rock. It is found that the pore fluid configuration within the individual pores and the pore geometries in which saturation conditions are varying are critical factors in determining the dielectric and elastic properties of partially saturated rocks. In addition, it is observed that the effect of saturation condition variations in a particular pore geometry increases as the pore shape becomes more crack-like. 11 The theoretical formulations describing the effects of pore-scale fluid distribution are used to analyze experimental data for a partially saturated tight gas sandstone and glass bead packing. Simple models that incorporate only the basic geometrical elements of the resulting pore-scale fluid distributions accurately predict the experimental data. It is also found that the same geometrical model can be used simultaneously to estimate the dielectric and elastic responses of a partially saturated porous medium. The effect of surface phenomena (e.g. electrical double layers and surface conduction) at the solid-fluid interface are incorporated into inclusion-based formulations of the dielectric response by employing the limiting case of a confocally-layered ellipsoid. It is found that the effect of surface phenomena varies as a function of both inclusion size and shape. in Table of Contents Abstract . ii Table of Contents iv List of Figures vii List of Tables xi Acknowledgements xii Chapter 1: Prologue 1 Chapter 2: A Sphere Packing Model Incorporating Contact Generation 5 Introduction 5 Undeformed Configuration 7 Overall Stress-Strain Relationship 8 Effective Moduli 15 Numerical Results 21 Conclusions 25 Chapter 3: Inclusion-Based Formulations for the Dielectric Properties of Heterogeneous Media 33 Introduction 33 Definition of an Effective Medium 34 Analysis of an Individual Inclusion 36 Classical Inclusion Formulations 40 iv Bruggemanâs Symmetric and Unsymmetric Fou1ations .44 Implicit Microstructural Relationships 48 Conclusions 51 Chapter 4: The Effects of Pore-Scale Fluid Distribution on Elastic and Electromagnetic Wave Propagation in Partially Saturated Porous Rock: A Theoretical Analysis 53 Introduction 53 Geometrical Description of Porous Rocks 54 Mathematical Formulations for Homogeneous Pore Spaces 56 Mathematical Formulations for Heterogeneous Pore Spaces 59 Numerical Simulations 65 Conclusions 71 Chapter 5: The Effects of Pore-Scale Fluid Distribution on Elastic and Electromagnetic Wave Propagation in Partially Saturated Porous Rock: An Analysis of Experimental Results 100 Introduction 100 Tight Gas Sandstone: Experimental Results 101 Tight Gas Sandstone: Numerical Simulation 105 Glass Bead Packing: Experimental Data 107 Glass Bead Packing: Numerical Simulation 110 Conclusions 111 Chapter 6: A Model for Incorporating the Effects of Surface Phenomena in the Dielectric Response of Heterogeneous Media V 131 Introduction .131 Response of a Homogeneous Inclusion with Surface Phenomena 132 Numerical Simulation 142 Conclusions 144 Chapter 7: Summary 155 References 159 Appendix 1: Constitutive Relationships and Generalized Permittivity 167 Appendix 2: Evaluation of the Depolarization Coefficient Associated with the Unique Principal Axis of the Inner Spheroid within a Confocally 171 Layered Inclusion vi List of Figures 2.1 Representative Sphere and Neighborhood 27 2.2 Point Contact Configuration 28 2.3 Near Contact Configuration 29 2.4 Computed Log (-e) as a Function of Log (p) 30 2.5 Computed Log (A.*) as a Function of Log (p) 31 2.6 Computed Log (*) as a Function of Log (p) 32 4.1 Computed 1c/1c(Sb=O.OO) for Homogeneous Pores 74 4.2 Computed c/G(Sb=O.0O) for Homogeneous Pores 75 4.3 Computed 1I1(Sb=O.O0) for Pores Filled with the Water-Wet Effective Pore Fluid 4.4 76 Computed a/(Sb=0.O0) for Pores Filled with the Water-Wet Effective Pore Fluid 77 4.5 Computed 1c/1c(Sb=0.00) for Water-Wet Confocally-Layered Pores 78 4.6 Computed a/a(Sb=0.00) for Water-Wet Confocally-Layered Pores 79 4.7 Computed 1c/1(Sb=0.00) for Pores Filled with the Air-Wet Effective 4.8 Pore Fluid 80 Computed 1c/1(Sb=0.00) for Air-Wet Confocally-Layered Pores 81 vii 4.9 Computed a/a(Sb=O.00) for Pores Filled with the Air-Wet Effective Pore Fluid 82 4.10 Computed a/a(Sb=0.00) for Air-Wet Confocally-Layered Pores 83 4.11 Computed k*/k*(Sbw=0.0O) for Homogeneous Pores 84 4.12 Computed 85 for Homogeneous Pores 4.13 Computed V/V(Sb=0.00) for Homogeneous Pores 86 4.14 Computed v:/v:(sbW=o.oo) for Homogeneous Pores 87 4.15 Computed k*/k*(Sbw=0.00) for Heterogeneous Pores 88 4.16 Computed Ii*II.L*(Sbw=0.00) for Heterogeneous Pores 89 4.17 Computed V/V(Sb=0.00) for Heterogeneous Pores 90 4.18 Computed Vâ/V(Sb=0.0O) for Heterogeneous Pores 91 4.19 Schematic Diagram for Varying Local SbW Level within the Pore Geometry Spectrum 92 4.20 Pore Geometry Spectrum for Model Sandstone 93 4.21 Computed 1c as a Function of SbW for the Model Sandstone 94 4.22 Computed 0 as a Function of Sb for the Model Sandstone 95 4.23 Computed k* as a Function of 5 bw for the Model Sandstone 96 4.24 Computed as a Function of Sb for the Model Sandstone 97 4.25 Computed V as a Function of Sbw for the Model Sandstone 98 4.26 Computed v as a Function of 5 bw for the Model Sandstone 99 5.1 Measured 1a of a Tight Gas Sandstone 115 5.2 Measured V, of a Tight Gas Sandstone 116 viii 5.3 Measured V of a Tight Gas Sdstone 5.4 Schematic Diagram of Pore-Scale Fluid Distribution during the . 118 Drainage Process 5.5 117 Schematic Diagram of Pore-Scale Fluid Distribution during the Imbibition Process 119 5.6 Pore Geometry Spectra for Models of the Tight Gas Sandstone 120 5.7 Computed ic for the Tight Gas Sandstone Model 121 5.8 Computed V for the Tight Gas Sandstone Model 122 5.9 Computed V for the Tight Gas Sandstone Model 123 5.10 Measured VP of a Glass Bead Packing MPa) 124 5.11 Measured V of a Glass Bead Packing (pe=lO. 33 MPa) 125 5.12 Measured V, of a Glass Bead Packing .Â°Â° 31 (Pe MPa) 126 5.13 Measured V of a Glass Bead Packing (pe=3 1.00 MPa) 127 33 . 10 (Pe 5.14 Pore Geometry Spectrum for Glass Bead Packing Model 128 5.15 Computed V for the Glass Bead Packing Model 129 5.16 Computed V for the Glass Bead Packing Model 130 6.1 Computed ic for models possessing a surface conduction (cz=1.000) 147 6.2 Computed ic for models possessing a surface conduction (a=0.100) 148 6.3 Computed ic for models possessing a surface conduction (=0.010) 149 6.4 Computed ic for models possessing a surface conduction (cz=0.001) 150 6.5 Computed a for models possessing a surface conduction (=1.000) 151 ix 6.6 Computed a for models possessing a surface conduction (cz=O. 100). 152 6.7 Computed a for models possessing a surface conduction Qx=0.01O) 153 6.8 Computed a for models possessing a surface conduction (x=0.001) 154 x List of Tables 2.1 Parameters for Sphere Packing Model 26 4.1 Dielectric Constants and Conductivities for Model Sandstone 73 4.2 Elastic Moduli and Densities for Model Sandstone 73 5.1 Local Level of Sbw Maintained in Pore Throats 113 5.2 Dielecthc Constant and Conductivities for Tight Gas Sandstone Model... 113 5.3 Elastic Moduli and Densities for Tight Gas Sandstone Model 114 5.4 Elastic Moduli and Densities for Glass Bead Packing Model 114 6.1 Dielectric Constant and Conductivities for Colloidal Suspension Model.. 146 xi Acknowledgements First and foremost, I wish to express my sincere gratitude to my research supervisor Rosemary Knight and my wife Zona for the support they have given to me over the last four and a half years. For this reason, I am deeply indebted to both of these very special individuals. I would like to thank the faculty, staff and graduate students of the Department of Geophysics and Astronomy, as well as the other members of the Rock Physics Research Program, for the many things (both tangible and intangible) they have graciously provided during my time at UBC. Further, I would like to acknowledge the contributions of Garry Clarke, Phil Loewen and Doug Oldenburg as members of my thesis supervisory committee. Without the guidance and counsel of three important individuals, my path in the field of geophysics, not to mention my life, would have been entirely different. In recognition, this thesis is dedicated to Lloyal Bacon, Pat Keiran and Mitch Yancey. Financial support was provided by the University of British Columbia in the form of a University Graduate Fellowship and an Izaak Walton Killam Predoctoral Fellowship. Funding for this research was obtained from the University Research Grants Program of Imperial Oil Limited. xii Chapter One Prologue âThe book of nature is written in the language of mathematics.â Galileo Galilei (1564-1642) Geophysics, as a discipline, is frequently concerned with the problem of determining the nature of porous rocks that cannot be directly observed. To achieve this goal, geophysical surveys are conducted to measure the spatial and/or temporal distribution of particular physical properties. However, these physical property measurements provide only an indirect means of determining and analyzing microscopic phenomena in porous rock. When one is concerned with these types of phenomena, it is necessary to establish an appropriate framework in which geophysical measurements can be interpreted. One approach is to perform laboratory measurements on porous rock samples that are believed to be similar to those that exist in situ. However, these types of measurements can prove to be of limited utility. Minor differences in terms of rock type and/or the conditions to which the rock is subjected can result in significant variations in physical behavior. Therefore, it is necessary to complement these laboratory results by developing mathematical formulations that describe the physical response of porous rocks in terms of its microscopic phenomena. In this study, theoretical formulations are developed which describe the role of certain microscopic phenomena in determining the observed physical 1 properties of porous rocks. Given the relative importance of electromagnetic and seismic techniques in geophysics, this study will be confined to the consideration of dielectric and elastic response of porous rocks. In unconsolidated sands, an important type of microscopic phenomenon is the interaction between neighboring sand grains. Experimental data for granular materials consistently violate the relationship between elastic wave velocities and confining pressure predicted by previous sphere packing models. Murphy (1982a) proposed that this discrepancy between these earlier models and the experimental results is due to the creation of grain contacts as confining pressure increases. A random sphere packing model that permits the generation of grain contacts is derived in Chapter Two. However, this type of formulation is specifically structured to incorporate interactions between neighboring particles. Hence, it is not well suited for the study of phenomena that are related to geometrical configuration of the constituents within a porous rock. In the range of frequencies of interest in the geophysical measurement of dielectric behavior (i.e. 100 kHz to 1 GHz), a dominant mechanism affecting the dielectric response of porous rocks is the charge polarization at interfaces on the microscopic scale (Poley et al., 1978). Further, it has been observed that the geometry of the pore spaces within a rock significantly affects its elastic wave velocities (Nur and Simmons, 1969; Toksoz et al., 1976). Therefore, it is necessary to use a formulation that incorporates sufficient information about the pore space geometries within the system to describe elastic and electromagnetic wave propagation in porous rocks. This is accomplished by employing an inclusion-based formulation to describe porous rocks. In this study, the porous rock system is viewed as a rock matrix into which inclusions representing the individual pore 2 spaces are embedded. In previous work, a diverse collection of mathematical expressions have been obtained when the inclusion-based formulation has been employed to describe a given physical property of composite materials. To clarify the relationship between these formulations, the major classes of the inclusion-based formulations are derived in Chapter Three in such a way as to avoid special cases that result from the use of specific inclusion types. The discussion is in terms of dielectric properties; however, the basic methods and principles used in the derivation are applicable to inclusion-based formulations for other physical properties, such as elastic wave propagation. A phenomenon that is characterized by the geometrical configuration of the constituents within a porous rock is the pore-scale fluid distribution within the partiallysaturated system. In Chapter Four, mathematical expressions are developed for the dielectric and elastic response of partially saturated porous rocks for various geometrical configurations of the pore fluids within the individual pore space. To determine the critical factors that control the relationship between pore-scale fluid distribution and macroscopic physical properties, a number of numerical simulations are performed utilizing these formulations. The basic results obtained from this theoretical study are applied in Chapter Five to the analysis of experimental data collected on partially saturated porous rocks. The measured apparent dielectric constant and elastic wave velocities for a tight gas sandstone undergoing an imbibition-drainage cycle (Knight and Nur, 1987a; Knight and Nolen Hoeksema, 1990) are modelled by means of simple geometrical scenarios for the resulting pore-scale fluid distribution. In addition, elastic wave velocity data for a 3 partially saturated glass bead packing obtained by Domenico (1977) are examined. In the derivation of the inclusion-based formulations, the existence of surface phenomena at the fluid-solid interface has been explicitly neglected, although such phenomena can have a significant effect on the electrical properties of porous rock. Previously, these phenomena were implicitly incorporated into mathematical formulations by assigning their effects to one of the constituents and adjusting its properties accordingly (Knight and Endres, 1990). However, this approach does not preserve the geometrical structure inherent in the system. In Chapter Six, an inclusion based formulation is presented that explicitly incorporates the effects of surface phenomena in the dielectric response of a heterogeneous medium. This is achieved by taking a limiting case of a confocally-layered ellipsoid as its outer shell becomes infinitesimally thin. 4 Chapter Two A Sphere Packing Model Incorporating Contact Generation Introduction The modelling of granular materials by means of sphere packings has been the subject of numerous papers. Many authors (Gassmann (1951a), Duffy and Mindlin (1957), Deresiewicz (1958), White (1965), and Walton (1975) ) employed regular packings (e.g. hexagonal close packing, cubic packing). Brandt (1955) presented the first formulation for the elastic behavior of random packings; however, only the overall volumetric changes and the effective bulk modulus of the medium were considered. More recently, Digby (1981) and Walton (1987) have presented more comprehensive treatments for the mechanical behavior of random sphere packings. In all of these studies, it was assumed that grain contacts were established in the initial, undeformed state and that no new grain contacts were generated during the compression of the medium. Further, some specific form of contact microphysics was always assumed. Experimental evidence has shown that these mathematical models do not adequately describe the observed elastic behavior of sphere packings. In all of the above models, the hydrostatic strain and the effective elastic moduli are shown to be dependent on the 2/3 and 1/3 power, respectively, of the hydrostatic confining pressure. This implies that the elastic wave velocities are related to confining pressure by essentially a 1/6 power law. Duffy and Mindlin (1957) presented experimental data for the elastic wave velocities of a 5 face-centered cubic packing of spheres. These data showed that the predicted elastic wave velocity was significantly higher than the measured velocity, and did not have a 1/6 power dependence on confining pressure. Further, as the allowable variation in sphere radii increased, the deviation between predicted and experimental velocity also increased. Domenico (1977) measured the elastic wave velocities of the random packing of air saturated glass beads under hydrostatic confining pressures between 2.7 and 34.4 MPa. A least squares fit to the data implied that both compressional and shear wave velocities have approximately a 1/4 power relationship with the confining pressure. Murphy (1982a) measured the compressional and shear wave velocities of a random packing of glass beads under vacuum and hydrostatic confining pressures between 0.1 and 35 MPa. A 1/4 power relationship was observed at low confining pressures, and a 1/6 power relationship was exhibited at high confining pressures. Murphy (1982a) stated that the 1/4 power relationship can be explained by the generation of contacts due to the closing of near-contact gaps and used a simple two-layer model to illustrate this point. However, a general formulation for a random packing was not derived. In this chapter, expressions are derived for the overall stress-strain relationship and the effective moduli of a granular medium under hydrostatic compression, which permit the generation of grain contacts due to the closure of near contact gaps. The mathematical formulation of Walton (1987) provides a basis which can be expanded to include the effect of contact generation caused by the closing of near contact gaps as the hydrostatic pressure increases. Further, Walton (1987) considers two cases of contact microphysics (i.e. perfectly smooth and infinitely rough spheres). The formulation presented here extends these results by allowing the use of general contact 6 microphysics. Undeformed Configuration Consider a random packing of identical solid spheres of radius R. Each sphere is elastically isotropic and homogeneous with LamĂ© constants .1 and . Let V be a volume within the packing containing N spheres and having dimensions which are very large compared to R. Within V, the random packing is assumed to be homogeneous. Let x (n) be the position of the center of the n th sphere relative to some given Cartesian coordinate system. The n th sphere is surrounded by neighboring spheres, some of which are in point contact with the n th sphere, while others are in near-contact, as in Figure (2.1). In Figure (2.2), the system of two spheres initially in point contact (in this case, the m th and n th spheres) is shown. Let the midpoint between sphere centers be denoted by 0, with a position vector x x (Â°) = - (Â°) given by (x (n) + x (m)) (2.1) . For this system, the inward unit vector along the line of the sphere centers, i is defined by I (nm) = (2R)-â (x (n) - x Cm)) . (2.2) The system of two spheres in near-contact is illustrated in Figure (2.3). In this case, let 0â 7 be the midpoint between the sphere centers, with a position vector x (Â°) = - (x (n) + x (Â°) given by (rn)). (2.3) Similarly, the inward unit vector along the line of the sphere centers, j for this system is defined by (nm) = [2(R+g(nm))]â (x (n) .. x (rn)), (2.4) where 2g (nm ) is the distance (or gap) between the surface of n th and mâ th spheres along the line of the sphere centers. Overall Stress-Strain Relationship In order to derive the overall stress-strain relationship, consider an initial deformed state, obtained by subjecting the boundary of volume V to a displacement u consistent with a uniform, compressive hydrostatic strain EI (where I is the identity tensor). Hence, for a position vector x in V, the displacement at this point if V were a continuum is given by n(x)=eI x=Ex . (2.5) It will be assumed that the sphere centers undergo a deformation which is consistent with this uniform strain and that all deformations are reversible. Then, 8 u (n) = ex (2.6) (n) defines the displacement of the n th sphere center. Similar expressions give the displacements for the centers of the m th and mâ th sphere. Consider the situation where two spheres are in point contact in the undeformed state. is the displacement of midpoint 0 with respect to the n th sphere center, then A(no) = 1 (u (m) - u (n)) = - RE i (nm) (2.7) As the hydrostatic strain increases, the contact area between spheres n and m becomes finite; and, a force, Fm), is exerted on the n th sphere due to its contact with the m th. Since the spheres are elastically isotropic and homogeneous, it can be assumed that Fm) is parallel to (nm) Hence, F () = i (nm), FN (2.8) where FN is the magnitude of the normal component of F(m) and &1m) defines the state of the contact between the n th and m th spheres. Since all original contacts undergo identical deformations for a given hydrostatic strain (nm) , = for all these contacts. For the spheres in near-contact, it is first necessary to establish a point contact before a force is exerted on the n th sphere. Since the relative displacement of one sphere center with respect to another is along the line of their centers, it can be shown that the point 9 contact is obtained when ii =g(m), ii (2.9) where A(fb) A (no =-(R+g(1m)) ei (nm), (2.10) ) being the displacement of the midpoint 0â with respect to the n th sphere center. Once this occurs, a finite contact area is established as the hydrostatic strain continues to increase; and, a force, F (nm) is exerted on the n th sphere due to its contact with the rn th sphere. Given the above assumptions, the force F(m ) is given by F (m) = FN = FN ( )) H (i i A (no) ()) H ( - (R + g (nm)) C I - - g (nm)) g (nm)) (nm) (nm), (2.11) where (nm ) defines the state of the contact between spheres n and mâ, and H is the Heaviside step function. Now that the forces on the n th sphere due to its contacts with its neighboring spheres are known, it is possible to relate the average stress within V to the average strain. The average stress, (ajj), within V is given by 10 (Gjj ) = J dV = â spheres where V is the volume of the n th sphere, J ydV, (2.12) V,i are the components of the Cauchy stress within the n th sphere, and the summation is over all spheres in V. For the n th sphere, the stress within V can be related to the tractions on S , the surface of the n th sphere, in the following manner: 1vGiJ J (x; dV + t)) ds, (2.13) = where is the traction on S and xâ = x x() is the position of a point on S with - respect to the center of the n th sphere. In general, it can be assumed that the contact areas will be small. This allows two simplifying assumptions. First, xâ may be approximated by xâ (x (m) - x (n)) = - R I () (2.14) for the original contacts and by R (x(m)xn)=Ri(â) xâ= 1 2 R+g() (2.15) for the generated contacts. Secondly, the integration of the i th component of the traction 11 over the contact surface is Fm) for an original contact and (hlm 1 F for the generated contacts. Using these approximations, Equation (2.13) becomes a(n) { ((nm) F(m + m) dV = x ) Ff )+ 1 F Fâ + ) }. (2.16) In this expression, the first term of the right hand side is the contribution of the original contacts and the second term represents the effect of created contacts. By using Equation (2.16) in Equation (2.12) and changing the summation from over the indices to over individual contacts, the following is obtained aj> =-. creat. { orig. cont. }. cont. (2.17) The fact that the medium is undergoing hydrostatic compression is now introduced. This means @ij)=Pijâ 12 (2.18) where p is the pressure applied to the surface of V and 5 is the Kronecker delta. By employing Equations (2.8), (2.11), and (2.18), Equation (2.17) becomes 3 i 6 P { = F(â )H + orig. cont (- (R + g( ))e - g(mn ))(i ))) creat. cont. }. (2.19) Since V contains a large number, N, of spheres, Equation (2.19) can be expressed in terms of average values for the type of contact involved as follows: p flc where (FN jj = (â)H â- { T FN () (j(nm) j(nm)) + (- (R + g(nm )) e - g(nm ))) ))) (2.20) and 11 are the average number of original contacts per sphere and the average number of near-contacts per sphere, respectively, and <...> represent average values taken over the contacts. Since it has been assumed that the packing is statistically homogeneous and isotropic on the scale of the dimensions of V, then original and created contact points have a uniform probability distribution over a sphere surface. In this 13 situation, Batchelor and OâBrien (1977) show that (iâ) j(mn)) = ) (j(nin (nm )) (2.21) = . Further, using =4R N 3 3(1 (2.22) , - where 0 is the porosity of the undefomied aggregate, Equation (2.20) becomes = 4R {lo FN () + (FN ((nm )) H (- (R + g(nm )) e - g(nm))) }. (2.23) The average value of the normal force at the created contacts can be determined if the probability density function of g (nm), f(g (nm)) is known. Since g (nm) , 0 and the Heaviside step function is zero if 4 g(nm)..Re(1+E) (2.24) , then, Equation (2.23) can be written as -RE(1 = 2 4R 1 +Ey FN() + FN() 14 f( g(mn ))dg(nmâ )}. (2.25) In order to obtain the overall stress-strain relationship, the form of the contact microphysics must be defined. A specific example is considered in the section on numerical results below. Effective Moduli In order to determine the effective moduli of the granular medium, it is necessary to consider the effect of the application of an additional incremental deformation to the sphere packing after the initial hydrostatic deformation. Therefore, assume that the boundary of V is subjected to a further deformation Ă¶uwhich is consistent with a uniform strain tensor Ă¶E. Then, the displacement of a point x in V, if it were in a continuum, would be given by u (x) = Ă¶Eâą x or , u (x) = & x. (2.26) It will be assumed that this incremental deformation is small enough that effects of contact generation and destruction can be neglected, that the sphere centers undergo an incremental deformation which is consistent with a uniform incremental strain, and that all incremental deformations are reversible. Hence, (n) = E . x (n) , or u = â. Given this as a starting point, then the incremental relative displacements of the midpoints 0 and 0â with respect to the n th sphere center are given by 15 (2.27) 6A(no) = - R Ă¶Eâą j(mu) (2.28) and - (R + g(nm)) 6E. j(mu)âą (2.29) Since the spheres are composed of homogeneous, isotropic material, it can be assumed that the incremental force at an original contact is given by = oF 0 F(âI) j (mu) + 0 F.(I) j (nm) (2.30) where (mu) nm) j(nm)) (rim)] j H -R[oE. j(nm)-(&Mimu)imu))j (rim)] II - â â R[OE . (mu) (231) - Similarly, the incremental force at a created contact is OF ={8FN(mu))j G) + 3FT(Imu))j (mu)} H((R+g(mu)) eg(mu)) where 16 , x (2.32) j (nm) = [SE. j(nm) II -(R+g(nm)) [SE. j(nm) -(R+ g(nm)) - - (s (& jnm)(nm))j (nm)] (233) (nm))j (nm)] 1 j(nm)i I and SFN and SFT are the magnitude of the normal and tangential component of the incremental contact force SF, respectively. The concept of the contact stiffness is now introduced. Let DN and DT be the normal and tangential contact stiffness, respectively. Then, DN=dâ dw (2.34) and DT=dT, (2.35) where w is the magnitude of the normal component of the relative displacement of 0 or 0â (i.e. parallel to i (â) or i (nm )) and r is the magnitude of the tangential component of the relative displacement of 0 or 0â (i.e. perpendicular to i (nm) or i (nm)). Since the incremental displacement is small, then SFN=DN Sw, and 17 (2.36) FT=DT &. (2.37) Hence, the contact forces can be rewritten as - { R (DN()-D T((I)) ()+D )(6 T((T)) Eâą i()) (2.38) for the original contacts, and m) =(R+g(nm)) D { T()oE (DN()-DT) () ()+ (e H(-(R+ g(mn)Eg(nm))) } (2.39) for the created contacts. Using the same procedure as was used to determine the relationship between pressure and hydrostatic strain, the following relationship is found between the average incremental stress, (&Yj), and the incremental average strain (ek1): = (i 0) 2 8itR DT () { lo ((&jk) R [ 2 (D N () (i() iâ)) + - DT (&) 18 ()) (&kli ShiI. ((t) (nm) j im)))] (nm). (nm). (nm) k 1 + 1 )+ [ )) 2 ((D N (iââ) $) i1 - DT ))) (R + g ( )) H ( ( i) + (DT(CInm)) ((&jk) (i1m) (R + g(mn - (R + g (rim )) H ( - (R + iâ9 + (&k) (i(â) ik(â)))] )) e - g(nm g (rim )) e - ))) (&) g(nm }. ))) x (2.40) Since the initial deformation was hydrostatic, the homogeneous and isotropic nature of the random packing is retained during the incremental deformation. Hence, Equation (2.21) and another result from Batchelor and OâBrien (1977), (rim) /.\ i 1 . (rim) j 1 . (rim) k 1 . (nm)\ 1 /.\ i 1 (mn) â / . (rim) j 1 . (rim k ) 1. (nm)\ / (2.41) may be used. By inserting Equations (2.21) and (2.41) into (2.40), the resulting expression can be compared with = câj where the effective moduli, Ck1 , (&M), can be expressed in terms of the effective LamĂ© 19 (2.42) constants, j.tâ and as follows = (2.43) ij 6 8 k1 +J.L (6ikĂ¶1 +6il6k). Hence, it can be found that = ((D N (m )) - DT { () (m (i() DT()) + - )) c - g (nm))) } )) H( - (R + g (mn ))e- g (mn ))) }. ))) (R + g ( )) H (- (R + g (nm , (2.44) and = ((D N (ci { (m )+DT (D N () + ))(R + g (mn D T (t)) + (2.45) Letting f(g (nm)) again be the probability density function of g (mn), then, as was done above, A. and can be written as 20 = -RE(1 +c) 1 j ( N { () (cifr )) - DT j(r (m N () ))) (R - DT ()) X + g (nm)) f (g (nm)) dg (nm) + (2.46) and = -Rc(1 +e)â j ( ((nm N 2 )) + { DT ! (D ((nmâ N () ))) (R + + DT ()) + x g (nm)) f (g (nm)) dg (m) } . (2.47) Numerical Results From the above analysis, it can be seen that the general elastic behavior of a contact generating granular medium is described by Equations (2.25), (2.46), and (2.47). In order to obtain numerical results, it is necessary to define the nature of the contact microphysics and f(g (nmâ)), the probability density function of g (nm) for the sphere packing. , The contact microphysics gives the form of the interaction between spheres in contact. Let us assume that spheres are infinitely rough (i.e. no slippage along the contact surfaces). Then, using the results of Walton (1978), the following expressions can be derived for the contact forces and contact stiffness: 21 (e)312 41 FN()= FN(cI )) = 4Râ 2 3tB 2 / 3 (- Re - g(nm )(i + e)) , (2.49) 2R(- c)1â2 tB (2.50) 4R(- c)â 12 iu(2B+C) (2.51) DN()= DT(â1)= (2.48) 3itB (2.52) iu(2B+C) (2.53) where B=âââ Ji 1 + (2.54) and c=âiâIâ 1 - (2.55) The nature of the probability density function, f(g (nm)) , has not been established. However, given the nature of the problem, it is plausible that g (mâ) has a Maxwell 22 distribution. Hence, f(g (nm)) is given by (Papoulis (1965)) 211(2 f(g(nm)) = 2 exp (g(Im)) 2 I 31 cc 2 (2 [ - (g(1m)) cc2) 4] H (g(11mâ)) (2.56) âą -1 The maximum of f(g (nm)) occurs at (Xâą 2lâ2 and has a value of 23/2 (cc e 1/2â).By using Equations (2.48) through (2.53), and (2.56), the following expressions can be * * derived for the confining pressure p and the effective moduli ? and II: ââ (i - - 312 E) B 2 3it J 21/2 RE 2(1+e) itâ 3 cc 1 io+ B(2B+c) 2 1o 21/2 } , (2.57) RE\ 1 x 3 lCcc31/21+81 { 1 r 1 (1 - 1 /2x2exp( x2)dx+( x) I * 11= (i - - E 1 +E) (1 - (1 - 0 o+1c - if 1 +E/ where 23 o , (2.58) , (2.59) } B(2B+C) 2 1O x)â /2x2exp( x2)dx+( x)â /2x3exp( x2)dx] ( )H2 0) (5B + C) (- 8 1 r I - x)3I2x2exp( x2)dx J (1Ă)C(c)h/2 *= Lo (1 21/2 itâ 3 cc 2 (_Re) x 3 1+e 1 (1 - x)â /2x3exp( x2)dx] } = (2)â (-RE 12 (2.60) . â1 +eI The integrals in Equations (2.57) (2.59) can be shown to reduce to generalized - hypergeometric functions (Gradshteyn and Ryzhilc (1965)). However, to evaluate these expressions, it is much more practical to perform a numerical integration. Table (2.1) gives the values of parameters used in evaluating Equations (2.57) - (2.59). These values are typical of those used for dense packings of glass beads in Domenico (1977) and Murphy (1982a). Values for r and 11 c were obtained from Bernal and Mason (1960). Figures (2.4), (2.5), and (2.6) illustrate log * * (- e), log ( ), and log (a ) as functions of log (p), respectively. In each case, a = 0.OO1R and 0.0001R were used. In addition, the results for the model of Walton (1987) are also shown for two cases: the number of contacts 1 = T + Tc and 11 = flo. In each case, the model of Walton (1987) gives linear relationships in terms of log (p), corresponding to the 2/3 power law for (- e) and the 1/3 power law for the effective moduli 2.*and From Figures (2.4) (2.6), it can be seen that three distinct regions - exist in the response of the contact generating model. At low pressures, the contact generating model asymptotically approaches Waltonâs model for 1 = o as p decreases. 11 At high pressures, the contact generating model asymptotically approaches the completely contacted model, i = r + câ as p increases. Between the two, a transition 11 zone exists. The onset of the transition zone as pressure increases is related to the mean value g (nmâ) (i.e. once a significant portion of the possible contacts are established, the transition begins.). This behavior of the transition is independent of the choice of the 24 probability density function, as it is a general response of the sphere packing due to contact generation. If an attempt were made to describe the response of the model in the transition zone by means of a power law relationship, Figures (2.5) and (2.6) show that the effective exponent would be greater than 1/3 for the relationship between the effective moduli and the confining pressure. Hence, the effective power law between confining pressure and elastic wave velocities would have an exponent greater than 1/6. This is consistent with the experimental results of Domenico (1977) and Murphy (1982a). Conclusions In this chapter, a model for the elastic behavior of a granular medium which allows the generation of grain contacts during hydrostatic compression is presented. The form of this model permits the use of a variety of contact microphysics and statistical distributions for the gap width of the near-contact. Numerical results for this model display significant differences in comparison with models that assume no contact generation, particularly in the transition region. Further, it can be observed that this transition zone will occur over the range of confining pressures of interest in geological and geophysical problems (106 i0 7 Pa) for values of a small in comparison to the - sphere radius. This means that extremely small near-contact gaps can have a significant effect on the elastic behavior of a granular medium. The behavior of the contact generating model is consistent with experimental results for random packings of glass beads. 25 A, = 2.1474 X 1010 Pa 2.9655 x 1010 Pa R = â 4.064 x i0 m 0.383 = 6.4 = 2.1 Table 2.1: Parameters used in numerical results. 26 / / \ / / / \ / I / I / \ \ / / / / / / V â â â â â Figure 2.1: Representative sphere in V (n th sphere), surrounded by neighboring spheres in contact (m th spheres) and near-contact (mâ th spheres). 27 Figure 2.2: The original point contact configuration. 28 Figure 2.3: The near-contact configuration. 29 â2.0 â15 C (U 0 â3.0 0 -J 3.5 / -4.0 / / A 4.5 5.0 6.0 5.5 6.5 7.0 7.5 Log(Pressure(Pa)) Figure 2.4: Computed log (-e) as a function of log (p) using the contact generating model where (1) a = 0.001 R and (2) a = 0.0001 R. Dashed lines are Walton (1987) model for rj rc, (Line A- A) and r = rio + lie (Line B Bâ). 30 - 8.2 A 8.0 7 / 7.8 0 -o 7.6 - E 1 C) 0 -J 7.4 7.2 / B, 7.0 A 4.5 5.0 6.0 5.5 6.5 7.0 7.5 Log(Pressure(Pa)) contact generating Figure 2.5: Computed log (A*) as a function of log (p) using the model where (1) cx = 0.001 R and (2) cx = 0.0001 R. Dashed lines are Walton (1987) model for ii = ro (Line A- Aâ) and ii Tb + Tic (Line B Bâ). 31 - 9.4 7 9.2 V V / V 9.0 / a D 8.8 - t3) 0 -j 8.6 7 / / 8.4 B7 A 4.5 5.0 6.0 5.5 6.5 7.0 7.5 Log(Pressure(Pa)) Figure 2.6: Computed log (j*) as a function of log (p) using the contact generating model where (1) a = 0.001 R and (2) a = 0.0001 R. Dashed lines are Walton (1987) model for 1 = lo (Line A- Aâ) and 11 =lo + Tc (Line B Bâ). 32 - Chapter Three Inclusion-Based Formulations for the Dielectric Properties of Heterogeneous Media Introduction In the range of frequencies of interest in the geophysical measurement of dielectric behavior (i.e. 100 kHz to 1 GHz), a dominant mechanism affecting the dielectric response of porous rocks is the charge polarization at interfaces on the microscopic scale (Poley et al., 1978). Inclusion-based formulations permit the incorporation of sufficient information about the geometrical configuration of the constituents in a porous rock to describe this phenomenon. Numerous studies have used inclusion-based formulations to investigate the relationship between microstructural features and macroscopic dielectric properties for heterogeneous media; these have been reviewed by van Beek (1967). All of these formulations describe the microstructure in terms of a background matrix into which inclusions are embedded; however, a diverse collection of mathematical expressions has been obtained. While some of this diversity can be attributed to the explicit use of specific inclusion types (i.e. inclusions with a given shape and composition), other variations result from fundamental differences in the methods used to derive these formulations. To illustrate these fundamental differences, as well as the common elements, the major classes of inclusion-based formulations are derived in this chapter. The degree of generalization that is maintained in this discussion is not found in previous work on the subject; therefore, the differences between the resulting expressions 33 are not caused by the explicit choice of specific inclusion types. In addition, this work has permitted us to identify a general class of iterative formulations and their associated hierarchical structure. It is important to note that the basic methods and principles used in the derivation of the dielectric properties are applicable to inclusion-based formulations for other physical properties. Definition of an Effective Medium To start the discussion, let us consider a heterogeneous medium that is composed of isotropic materials possessing causal, linear electrical rheologies that are described in terms of generalized permittivities (see Appendix 1). Let a volume V be some representative volume element within this heterogeneous system which has dimensions that are large in comparison to the scale of the inhomogeneities, yet small in comparison to the size of the total system. The material in V is subjected to an applied uniform electric field that is time harmonic with angular frequency . It will be assumed that the wavelengths and attenuation lengths of this field are large in comparison with the dimensions of V. Therefore, a quasi-static approach is appropriate for determining the propagation of electromagnetic waves within this volume. The microscopic structure is uniformly distributed such that the system is homogeneous and isotropic on a macroscopic scale. Hence, the macroscopic electrical behavior of the heterogeneous system can be simulated by the response of an effective homogeneous material with a generalized permittivity . Let D ard E be the generalized electric displacement and electric field strength within the system, respectively, in complex notation. Then, the effective generalized permittivity of the 34 material in V is defined by (3.1) where D(x)dv (32) E(x)dv (3.3) Jv and 1 E=V f iv are the volume averages of D ani E respectively. In addition, the pointwise , constitutive relationship (3.4) D(x) =eg(x)E(x) is valid at all points x in V. Rearranging terms in Equation (3.4) and taking the volume average, one obtains + v.âJ [cg(x) 35 - E] E(x) dv. (3.5) Given the definition of the effective generalized permittivity (Equation 3.1), Equation (3.5) implies that vâ [eg(x)-egjE(x)dv=O. (3.6) iv Hence, the volume average of the pointwise difference in polarization between the heterogeneous and effective media must vanish. This equation defines the fundamental relationship between a heterogeneous system and its equivalent effective medium. To evaluate Equation (3.6), the structure of the underlying medium used in inclusion-based formulations will now be considered; it will be seen that it is necessary to analyze the response of the individual inclusions in order to accomplish this procedure. Analysis of an Individual Inclusion It will be assumed that the microstructure within the heterogeneous medium can be described in terms of a homogeneous matrix with a generalized permittivity Eg(0) into which homogeneous inclusions composed of materials having generalized permittivities are embedded. Considering the structure of this system, Equation (3.6) may be rewritten as V ; VegÂ°+ J (eg0 eg(Â°))EO)(x)dv vinc 36 (3.7) where is the volume of the inclusions in V and E ( (x) is the electric field strength 0 within the inclusions. The two terms on the right hand side of Equation (3.7) can be viewed as a reference and perturbation term, respectively. In general, it is very difficult to precisely determine the electric field within the inclusions. However, it is possible to estimate E (x) within an individual inclusion by treating it as a single inclusion embedded in an infinite, homogeneous background with a generalized permittivity egO). This system is subjected to an applied field with uniform strength Ea at infinity. The values of EgÂ° and Ea are selected such that interactions with the other inclusions in the system are simulated. Given the macroscopic isotropy of the heterogeneous system, it will be assumed that Ea and E are oriented in the same direction. The quasi-static behavior of this system is govern by cgO) V i.i,(x) 2 0 (3.8) 0 (3.9) within the background matrix and )v 4 ,() (x) 4 2 inside the inclusion. The term qi = is the electric potential; and, E (x) = V - i.ji (x). Continuity of the electric potential and the normal component of the generalized electric displacement lead to the following boundary conditions along the surface S that forms the interface between the inclusion and the background: 37 jj(b)(x) =v()(x) Is 5 I (3.10) and ([egOâ) =([eg(O VNI(â)(x)] ii) Is, 5 ii) 1 r0â)(x)] 1 v (3.11) . where ii is the unit normal vector on S directed into the background. Such boundary conditions exclude the existence of surface phenomena, such as an electrical double layer, along the surface S. (x), the inclusion shapes must be 0 In order to obtain an analytic solution for E ( restricted to ellipsoids. The ellipsoid under consideration has semi-axes of lengths , and a3, respectively, where a 2 1 ai, a a2 . The shape of the ellipsoid is uniquely 3 a defined by the aspect ratio pair (xl,c2) where , (3.12) Using the results of Stratton (1941), it can be shown that the electric field within this ellipsoid is uniform and has a strength given by )] cos 2 A(xi,cL E ()(x) = Ea8g , {[egb +(egO) egOâ)) 4 0 - where Ea is the magnitude of Ea, 0 j is the angle between the j th semi-axis and Ea, ) is the depolarization 2 is the unit vector along the j th semi-axis. The term A(cq,c 38 (3.13) and coefficient associated with the j th semi-axis and is uniquely determined by the shape of the ellipsoid. In terms of the semi-axis lengths, the depolarization coefficients are defined as rÂ° A(al,cx2)=-ala2a3 J r [ 4 (s+a2) J] i-i (2)j (j=1,2,3). ds 2 (3.14) k=1 Jo The depolarization coefficients for the three principal axes are related by ( A ( ) 2 ( 3 , 1 ) +A cx x = x a 1. (3.15) In order to preserve macroscopic isotropy in the heterogeneous medium, identical ellipsoidal inclusions are assumed to be oriented in a uniformly random manner. The electric field strength inside such an inclusion is estimated by averaging the result given by Equation (3.13) over all possible orientations of the inclusion with respect to Ea. Performing this process, it is found that the expected value of the estimated electric field strength within an inclusion is (b) (E (Ox) { 3 [eg(b)+(eg()_eg(b) )] 2 , 1 A( } Ea (3.16) This quantity will be used in Equation (3.7) to approximate the electric field inside the inclusions within the heterogeneous medium. It will now be shown that major classes of inclusion-based formulations can by defined by the manner in which 39 cgOâ) and Ea are used to specify the inclusion interactions. Further, no explicit restrictions will be placed on the shape and composition of the ellipsoidal inclusions permitted in the medium; thus, the expressions which are obtained are of a more general nature than those previously given in the literature. Classical Inclusion Formulations Let us consider a heterogeneous medium in which a continuous distribution of ellipsoidal inclusion types are present. It will be assumed that the nature of the inclusion interactions is such that the same values of cgOâ) and E a can be used in the determination of E(x) within the individual inclusions. Using the result given by Equation (3.16), Equation (3.7) becomes e;E=c10)E+1e0)EaJ Egâ JJ 1 a1,c2) cjfl(4 , 0 x Z2 1 ,c12) th 2 f{egO, cg(Â°), eg(O, a 1 deg(â), da (3.17) where 1 ,a2)] f(egO), g(O), Eg, 0i ,a2) =(eg(â) eg(O)) [Eg +(g(â) egOâ)) A (a - - and E is the magnitude of L The distribution of inclusion types in V is described in 40 (3.18) terms of the fractional volume density function c(e),x1,cz2). For a distribution which is discrete in any inclusion parameter (i.e. cg(i), (X , (X2), 1 the appropriate integral is transformed into a summation and the definition of cj is modified. Equation (3.17) represents a group of formulations which will be referred to as classical inclusion formulations; various subdivisions of this group are defined by the values chosen for eg(b) and E a in order to simulate inclusion interactions. The magnitude of these interactions is dependent on the inclusion types present and their concentrations, as well as the relative contrast in the generalized permittivity of the constituents. If the inclusion concentration is very dilute, reasonable accuracy is obtained by assuming that there are no inclusion interactions (i.e. eg(1) eĂ§o) and E a = E) (Wagner, 1914; Burger, 1919; Fricke, 1953). Hence, Equation (3.17) can be expressed as = [1 + - J ff tZ f(eg(Â°),EgÂ°),eg,a1,a2) , (x2) 1 cinc (eâ, a x Z2 2 da da 1 dcg(1)] , (3.19) which is a generalization of what has previously been referred to as the dilute approximation. At higher levels of inclusion concentration, it is necessary to incorporate inclusion interactions into the analysis. One approach for simulating this effect is to assume (b) (o) and determine how the value of E a is modified by increasing inclusion concentration. One method, which has its origins 41 in the work of Lorentz and others (Mossotti, 1850; Clausius, 1879; Lorentz, 1880; Lorenz, 1880), achieves this by equating the sum of the induced dipole moments due to the inclusions in a given spherical volume of the heterogeneous medium with the induced dipole moment of an equivalent sphere composed of the effective medium (Landauer, 1978). From this analysis, it can be shown that *+2(O) g Ea E. (3.20) 3E.âLÂ°i Hence, equation (3.17) can be expressed as E*E(0) (0) +8 2 g f (() = 9 ( I ( ( I I JE Jai cj (e),al,a2) x J2 eg(O), CgXl X2) 2 do dx 1 dCg, (3.21) which is a generalization of what is commonly referred to as the Lorentz approximation. In addition, several other techniques for estimating the apparent field intensity E a have been suggested. Fricke (1924,1953) presented an approach where E a is assumed to be equal to the average electric field strength in the matrix material. Sihvola and Kong (1988) gave a technique where E a acting on a given inclusion is the sum of E and the strength of the depolarization fields due to the other inclusions. An alternate approach for simulating inclusion interactions is to assume that E a = E 42 and determine the appropriate value of egOâ) that describes the effective background experienced by the inclusions. From an analysis of experimental results, de Loor (1953, 1964, 1968) determined that egOâ) has a value between e and Â°) in lossless composites and that no unique value of egcâ) adequately describes all heterogeneous media. Of (BĂ¶ttcher, 1945; Polder and van Santan, particular interest is the case where egOâ) = e 1946; Taylor, 1965), which is commonly referred to as the self-consistency condition. This value of egOâ) was determined by BĂ¶ttcher (1945) on the basis of the theory of reaction fields (Onsager, 1936; BĂ¶ttcher, 1942). Incorporating this condition into Equation (3.17), one obtains = (o) + J Jâ Eg aj 11 (e cj )c1c2) 0 X âą. X2 f(e,egÂ°,eg,cx1,x2 1 degâ, 2 da ) dO (3.22) which is a generalization of the average field or self-consistent approximation. All classical inclusion formulations are obtained by determining the perturbation term in Equation (3.7) with respect to the background matrix e0). Hence, the accuracy of this group of formulations can be limited in terms of inclusion concentration. For cases of high levels of inclusion concentration, another group of formulations based on the two methods proposed by Bruggeman (1935) is commonly used. 43 Bruggemanâs Symmetric and Unsymmetric Formulations Bruggemanâs unsymmetric approach, commonly referred to as the differential effective medium or iterated dilute approximation (Norris et al., 1985a; Avellaneda, 1987), utilizes an iterative construction process to fabricate the effective medium. Starting with a homogeneous background matrix with generalized permittivity (0), each iteration involves the embedding of an infinitesimal portion of the inclusions present in the heterogeneous system. After each embedding, the dilute approximation (Equation 3.19) is used to determine the effective properties of the resulting medium; this resulting medium is used as the effective background matrix in the next embedding step. Using the variable volume process of Norris et al. (1985a), the generalized form of 0 be the initial volume the iterated dilute approximation can be derived as follows. Let V of the system prior to the homogenization process; then, the volume of the system after the homogenization process (Vf) is given by Vf= i v ( 0 )-1 - (3.23) where Ăžinc = J J jf 0 âŹgâ 1 2d 1 deg(â) dt ) dx 2 cj (t;e(),cti,cx (3.24) 2 is the total fractional volume of the inclusions. The fractional volume density function 1 cj (t;e), cx , CL2) is expressed in term of a homogenization parameter t that specifies the 44 order in which inclusions are incorporated into the effective medium during the homogenization process. For brevity, the term within the brackets in Equation (3.24) will Therefore, the total volume of the heterogeneous medium at any be denoted by point in the homogenization process, as specified by the value of t, is given by 1 V(t) = (i - 0 )â V 1 (s) ds (3.25) . Let us now consider the change in the generalized effective permittivity of the system due to the embedding of a very small volume of inclusions which denotes the increment from t to t + At in the homogenization parameter. If e(t) and E(t + At) denote the general permittivities of the effective background matrix and resulting medium at this embedding step, then the dilute approximation (Equation 3.19) gives v(t+At) e(t+At)=V(t+At) E;()+iE;()>< { fI (lĂ)-1J 0 v cxj t 2 )da f(c(s),e(s),4),xi,c 2 do 1 ,xi,a 0 ci(s;c ) 2 a2 deg(1)] ds } . (3.26) As the volume of the embedded inclusions becomes infinitesimally small (i.e. At â*0), then the change in c due to the homogeneous process is described by the differential 45 equation de(t) - .11 .flds]e(! =[l J JJ 1 a c2 ),xi,a 0 (t;c ) f(e (t),c (t),e,al,x2 1 d4â) 2 dc )dc . (3.27) 2 a The generalized permittivity of the heterogeneous medium is determined by integrating Equation (3.27) from t =0 to t = 1 with the initial condition e(t =0) = eg(Â°). Variations in the embedding sequence can result in significantly different values of e (Norris et al, 1985a). This initially appears to lead to ambiguities in the application of the iterative dilute approximation. However, it will be shown later that the nature of the microsiructure within the heterogeneous medium uniquely defines this homogenization process. In addition, it should be noted that the embedding sequence determines the nature of effective background matrix for a particular inclusion as the value of vanes between * . Eg(â) . âand Cg dunng the homogenization process. Bruggemanâs symmetric effective medium theory, which is also referred to as the coherent potential or effective medium approximation (Landauer, 1952; Stroud, 1975; Milton, 1985), carries inclusion-based formulations to their logical extreme by considering the heterogeneous medium to be an aggregate of inclusions. This structure implies that Equation (3.6) can be rewritten as 46 - ejE ()(x) dv =0. (3.28) Assuming E a = E and the self-consistency condition, one obtains J JJ ai )f(E, 2 Cinc (cgâxi, cX 1 deg() = 0 e, e,cxi,cx2 )dcx 2 da (3.29) a2 where I II (3.30) ),cxi,a2)da2dxidEg cinc(eg 1 0 . 2 a Norris et al. (1985 a, b) and Avellaneda (1987) have established relationships between Bruggemanâs symmetric and unsymmetric approaches. In particular, these results imply that the value of e, obtained from the iterated dilute approximation as Ăj case when all embedding steps are identical (i.e. c 1 converges to the Eg* (t, , 1 cx cx2) = 1 c â 1 for the (g(), , 1 cx )) 2 cc given by the coherent potential approximation. Further, Milton (1984, 1985) has shown that the coherent potential approximation is physically realized by a homogenization-type process. This suggests that a generalized version of the effective medium approximation can be obtained by permitting cj to vary during the embedding process while maintaining Ăj= 1. At this point, the major classes of inclusion-based formulations have been derived in their generalized form. The question that now arises is which of these expressions is the 47 most appropriate for a given heterogeneous medium. One criterion is to compare the microstructure of the heterogeneous medium with the underlying microstructural model for a given class of theoretical formulations. Implicit Microstructural Relationships Each theoretical formulation has an underlying microstructural model that determines the applicability of the formulation to a given heterogeneous medium. While some features of the microstructure are explicitly specified (e.g. designation of a background phase and specification of inclusion types), other relationships are implicit in nature. Of particular interest are the implicit microstructural relationships which result from use of the Bruggeman techniques for higher levels of inclusion concentration. The formulations derived in the preceding sections exhibit no explicit dependence of the dielectric response on inclusion size. In fact, it has been proven that the electrical properties are invariant with respect to the size of the inhomogeneities in the system under the conditions assumed in this study (Cohen, 1981). However, an analysis of experimental data for a number of simple systems indicates that there is a relationship between the distribution of the inclusion sizes present within the system and the accuracy of a given theoretical formulation (Pearce, 1955; De La Rue and Tobias, 1959; Meredith and Tobias, 196 1,1962). It can be argued heuristically that individual inclusions of a given size are embedded in an effective background material composed of the original background matrix material and the smaller size inclusions. Therefore, an iterative construction process will be used 48 to determine the effective generalized permittivity of the system. Starting with an original background matrix inclusions are sequentially embedded in order of increasing size. At the k th iteration, all inclusions of a given size are embedded and the effective generalized permittivity of the resultant system is determined; this system becomes the effective background matrix for the (k+1) th embedding step (i.e. c (k) = C(b)(k+ 1)). Let us initially consider a system containing K discrete sizes of inclusions. The variable volume process (Norris et al., 1985a) will be used during the iterative embedding procedure. Using a classical inclusion-based formulation, the embedding of the k th size inclusions is described by + a j J JJ 2 a cj (k;eg(), i, k+1 f JJ )da (j;c(),ai,a d ai de(1)] cj 2 2 a ) f(e(b)(k)e(b)(k) 2 da 1 dEgâ, Eg(â)ai cL2) da (3.31) where cj (),a a 5 e , ) is the volume fraction distribution function of inclusions of the 2 (k; 1 11 k th size and E a is determined by the approximation used. The effective generalized permittivity of the heterogeneous medium is given by e = Let us consider the case when a piecewise continuous spectrum of inclusion sizes is present in the system. This implies that the construction process involves the incorporation of an infinitesimally small inclusion volume during each embedding. In this case, the construction process described by Equation (3.31) becomes 49 dE(s) = eâ(s) J Jf 1 a where smin and max 5 { JSmax [J .L I dX 2 1 dEg(i)] dt} cj (t; C),Xi,cX2) dx >< cj (s;c,a1,x2)f(e(s),E(s),eg ,a1,cz2 )dcc2da1deg 0 (3.32) 2 a are the minimum and maximum inclusion sizes present in the medium. The effective generalized permittivity of the heterogeneous medium is obtained by integrating this initial value problem from E(s)=Eg(0) to E(5max) = Comparison of this result with Equation (3.27) shows the equivalence between this type of hierarchical medium and the iterated dilute approximation. Further, this equivalence demonstrates that the homogenization process in the iterated dilute approximation is uniquely defined in terms of the microsiructure of the medium (i.e. the embedding sequence of the inclusions proceeds in order of progressively larger size). However, the heuristic argument used to motivate the hierarchical models described by Equation (3.31) implicitly assumes that the inclusion of a given size are significantly larger than those incorporated in the effective background matrix. This implies that the relative concentration of inclusions having similar sizes is sufficiently dilute such that their mutual interactions can be neglected. For high levels of inclusion concentration, this implies that a very broad range of sizes must be present. This conclusion is supported by experimental observation (De La Rue and Tobias, 1959). The nature of the self-consistency condition used in the average field and coherent potential approximations now becomes apparent. For all inclusions to experience the 50 effective medium e as a background matrix, the distribution of inclusion types must be ),x1,c2)). This implies that the 0 identical for all sizes (ie.cinc (s; Eg,x1,cx2)=cmc (cg inclusions must be present in a wide range of sizes and fill the entire space (ie.pjn i). The generalized version suggested for the effective medium approximation differs from the coherent potential approximation in that distribution of inclusion types varies with size; hence, this generalization does not utilize the self-consistency condition. Another implicit relationship of the Bruggeman formulations is the degree of connectivity a given constituent possesses in the underlying model. To investigate this relationship for the coherent potential and iterated dilute approximations, Yonezawa and Cohen (1983) and Norris et al. (1985 a,b) employed simple two-component systems consisting of conducting and insulating phases to determine the percolation threshold of the resulting medium. In the iterated dilute approximation, the designated background component remains connected over all levels of inclusion concentration for nondegenerate cases of inclusion shape (i.e. needles and flat disc). In contrast, the coherent potential approximation always exhibits non-zero percolation thresholds beyond which complete connectivity of a given component ceases to exist; the value of this threshold is a function of the inclusion shapes. Conclusions The generalized versions of the major classes of inclusion-based formulations used to determine the macroscopic physical properties of heterogeneous media have been derived in this chapter. All formulations of this type originate from the requirement that the pointwise deviation in behavior between a heterogeneous system and its equivalent 51 effective medium must vanish in some average sense. To evaluate this criterion, it is necessary to determine the responses of the individual inclusions; their behavior is estimated by analyzing single inclusions embedded in an infinite, homogeneous background. Inclusion interactions are simulated by means of the imposed conditions used in this analysis; the major classes of inclusion-based formulations can be defined on the basis of the technique used to simulate these interactions. Each class of inclusionbased formulations has an associated microstructure that specifies the topological relationships between the various components. Of particular significance are the hierarchical structures, in terms of inclusion sizes, that are implied by use of an iterated embedding process. Further, the physical realization of formulations employing the selfconsistency condition represent a specific case of these hierarchical structures. In the next two chapters, inclusion-based formulations will be used to investigate the effects of variations in microstructure due to changes in pore-scale fluid distribution on the dielectric and elastic responses of partially-saturated porous rocks. 52 Chapter Four The Effects of Pore-Scale Fluid Distribution on Elastic and Electromagnetic Wave Propagation in Partially Saturated Porous Rocks: A Theoretical Analysis Introduction It is well established that the elastic and electrical properties of rocks vary as a function of the level of water saturation (SW) in the pore space. Hence, it is assumed that the level of water saturation within a porous rock can be determined by measuring these properties. However, it is now apparent that the specific form of the functional relationship between a particular physical property and S is significantly affected by the pore-scale fluid distribution occurring within the medium. This has been observed in laboratory data for elastic wave velocities (Domenico, 1976, 1977; Knight and Nolen Hoeksema, 1990), elastic wave attenuation (Bourbie and Zinszner, 1984), electrical resistivity (Longeron et al., 1989), and dielectric constant (Knight and Nur, 1987a). Significantly different values of a given physical property were observed for the same level of water saturation in these studies by varying the pore-scale fluid distribution through the use of differing saturation methods. Many of the theoretical formulations commonly used to model the physical properties of partially saturated rocks cannot account for this variation in functional dependence with changes in pore-scale fluid distribution. This is due to the inadequate description of the geometrical configuration of the constituents within the medium. In this chapter, 53 mathematical formulations are derived that allow the incorporation of the necessary geometrical information about pore-scale fluid distribution. These formulations are used to undertake a theoretical analysis of the effects of pore-scale fluid distribution on elastic and electromagnetic wave propagation in partially saturated porous rocks. Geometrical Description of Porous Rocks Let us consider the description of the geometrical configuration of the constituents within a porous rock. In order to utilize the inclusion-based formulations discussed in Chapter Three, this system will be viewed as being composed of a rock matrix into which inclusions representing individual pore spaces are embedded. By using this discretization of the system, it is possible to specify the saturation conditions within the individual pore spaces; however, it also implies that intra-pore phenomena control the physical properties of interest. For elastic and electromagnetic wave propagation, communication between individual pores is minimized when sampling frequencies are sufficiently high and/or the permeability of the pore system is adequately low. Under these conditions, pore connectivity can be neglected in the microstructural description of a porous rock when determining its dielectric and elastic properties. Since the volume fraction of the pore space is usually small, accurate results are commonly obtained by using the classical inclusion formulations with this form of discretization. In order to avoid possible concentration limitation associated with the dilute approximation, the Lorentz approximation will be used to compute the dielectric and elastic response of the porous rock system. The differences between results obtained from use of the Lorentz approximation and the higher order iterative dilute 54 approximation are small for parameters that realistically describe porous rocks. Information from pore cast (Swanson, 1979) and SEM observations (Timur et al., 1971; Brace et al., 1972; Sprunt and Brace, 1974) indicates that the pore space geometries can be adequately described in terms of oblate spheroids. The shape of an oblate spheroid is uniquely determined by a single aspect ratio a (i.e. the ratio of the minor semi-axis length to the major semi-axis length) that has a value between 0 and 1, inclusively. The description of the pore space geometries used in this analysis will be restricted to these shapes. Hence, the pore shapes used in this description range from spheres representing the major equidimensional pore volumes to discs representing cracks. By adopting this restriction on the permissible pore shapes, the computational effort required to determine the physical properties is significantly reduced. In particular, the depolarization coefficients which appear in the inclusion-based formulations for the dielectric response can be expressed in a simple analytical form. For a non-degenerate oblate spheroid (i.e. 0<a<1), the depolarization coefficients are given by (Landau and Lifshitz, 1960) (e...tan-le) 2 (=A(a)= 3 A ) 2 ai,cc 1Â±e (4.1) A ) ( , 1 ) 2 =,-(1 =A a a -A(a)) (4.2) and 55 where e = (y2 - 1)112 (43) For the degenerate cases of a sphere (i.e. ct=i) and a flat disc (i.e. cic=O), the depolarization coefficient associated with the unique principal axis is and 1, respectively. The shapes of the individual pore spaces are assumed to be independent of changes in the local saturation conditions. Hence, the distribution of pore geometries within a given rock is considered to be invariant and will be described in terms of a pore geometry spectrum cpore(x). While cpore(ct) can be viewed as a continuous spectrum, a discrete spectrum consisting of M differently shaped oblate spheroids will be used to describe the pore space geometries. Using the above description of the porous rock system, the mathematical expressions for the dielectric and elastic response of partially saturated porous rocks will now be determined. The forms of these expressions depend on the nature of the pore fluid configuration within the individual pore spaces. While porous rock can contain numerous pore fluids, the systems considered in this study will be assumed to contain only two distinct pore fluids. These systems are sufficient to illustrate the influence of pore-scale fluid distribution on the macroscopic physical properties. Mathematical Formulations for Homogeneous Pore Spaces Let us initially consider a situation where the incliviĂ§lual pore spaces are completely 56 filled with one of the two pore fluids; this implies that they are represented by homogeneous inclusions. Using Equation (3.21), the Lorentz approximation for the dielectric response of this system is given by * (0) gg( ) 0 5 C 2 Eg+ r2 1 M (4.4) pore(Xm) n=1 = where g(0),4fl)m) = (n(o)) { [e0)+(E(n)-e(0)) A(Xm)] â+ 4(e(Â°)+$ââ)+(E(Â°)E)) A(Gm)] In this expression, 0) -1 } âą is the generalized permittivity of rock matrix; and, (4.5) 4 (n= 1, 2) are the generalized permittivity of the pore fluids. Further, the saturation level of the n th pore fluid in the m th shaped pores is given by where s)+s)=i (m=1,...,M). (4.6) Analogous expressions for the effective elastic moduli of this system can be obtained from the formulation derived by Kuster and ToksĂ¶z (1974). The effective bulk modulus 57 k* and effective shear modulus k* -k(Â°) 3k*+40) â â 1 2 -[n=1 are given by M m=i s)cre(cxm) ( k)-k(Â°)(O))TUJJ (k(Â°), (Â°), k(n), (n), cLm)] (4.7) k(o) and 1 f 2 M sc(c) x n=1 m=1 FT â- 3 T .](lc0),.t(0),1 ),pâ),Xm) } , (4.8) respectively. The elastic moduli of the rock matrix are denoted by k(Â°), (o) and those of the pore fluids are represented by k(h1),.L(n1) (n=1, 2). The functions Tjj j and[T. 3 ijij i 3 are complicated expressions involving x and the constituent elastic moduli; the reader is referred to Kuster and Toksoz (1974) where both functions are defined. The effective density p* is given by 2 ,* M (i-ppoe)p(0)+ s)cpore(xm)p(â) (4.9) n=1 m=1 where M (Ppore cpore(czm) m=1 58 (4.10) is the total porosity in volume fraction, p(O) is the density of the rock matrix, and pOt) (n=1, 2) are the densities of the pore fluids. Once the effective elastic moduli and density are determined, the effective compressional and shear wave velocities are obtained from 1 (4.11) p and v: = L 2 (4.12) p respectively. To this point, only pore spaces that are completely filled with a single pore fluid have been considered. However, it is probable that several distinct pore fluid species will simultaneously occupy the same pore space within a partially saturated porous rock. Let us now consider the dielectric and elastic response of a porous rock when this situation occurs. Mathematical Formulations for Heterogeneous Pore Spaces When several immiscible pore fluids coexist within an individual pore space, it is necessary to use a heterogeneous inclusion, i.e. an inclusion containing two or more constituents, in order to describe the configuration of the pore fluids. In the case of the 59 elastic properties, the contents of a heterogeneous pore space can be replaced in the theoretical formulations by a homogeneous inclusion filled with an effective pore fluid regardless of the geometrical configuration of the fluids within that pore space. To show this, let us consider an individual pore space with a total volume V occupied with volumes V 1 and V 2 of the inviscid fluid species 1 and 2, respectively. Only the hydrostatic pressure component of the stress tensor exists in either of the fluids. Assuming that surface tension can be neglected, continuity of tractions at an interface between the fluids implies that P1 = P2 (the pressure in the respective pore fluid species). Let ap be an incremental change in the hydrostatic pressure of the pore fluids, then 1 V (4.13) and (4.14) where V 2 are the incremental change in the volumes occupied by fluids 1 and 2, 1 arxl V respectively. The contents of the pore volume are replaced with an effective fluid with bulk modulus kc and shear modulus = k(e) =0, then = 1 + aV (av ). 2 Rearranging the terms in Equation (4.15) and using Equations (4.13) and (4.14), we 60 (4.15) obtain U))) 5 ( (k(e) where (1) 1 = s(â) (k(â) ) â + (1_s(1) ((2) ) ) â (4.16) , is the volume fraction of pore fluid 1 in the pore space. It should be noted that given the above assumptions, Equation (4.16) is valid for any configuration of pore fluids within the pore space. To determine the elastic properties of the medium, the heterogeneous inclusions are replaced by homogeneous ones filled with the appropriate effective pore fluid. While a continuous spectrum of (i) values can occur within the heterogeneous pores of a given shape, it will be assumed in this study that all heterogeneous pores of a given shape have the same saturation level )âą Hence, the effective elastic moduli are given by k* k(Â°) 3 k*+ 4 â () 1 F k(e)(5(1)) Cpore (cCm) â 3k(Â°) 1 (k.(Â°), ,(Â°), (e)($)), (e)($)), am) Tjj - - k(o) 4(O) ] (4.17) and * 6 M (0) (k(Â°) +2 (o)) + (Â°) (9 k(Â°) + 8 (o)) 61 = cpore (am) x I(O)(3k)Â± 4i.(O)) 3 [T - - ijj] (k(O), p(O),k(e)(5j)), (e)(S)), (Lm) 1 T } . (4.18) In contrast to the elastic behavior, it is necessary to specify the geometrical configuration of the pore fluids within the heterogeneous pore spaces when mathematical expressions for electromagnetic wave propagation are derived. Two configurations of pore fluids for which analytical descriptions of the response can be obtained will now be considered. In both cases, 40 aal 42) will denote the generalized permittivity of the non- wetting and wetting pore fluids, respectively. First, let us consider the case where an effective pore fluid is generated by embedding small spherical inclusions of the nonwetting phase into a background of the wetting fluid such that the resulting assemblage resembles a gas bubble-water mixture. Using the Lorentz approximation, the generalized permittivity of the effective pore fluid 4e) is given by C(e) ((i) g - E(e)(S(1)) + (2) g = 2e2) () (i) g - (2) g , + (4.19) 2Ev) where s(1) is the fractional volume of the non-wetting phase in the pore space. Equation (4.19) is valid for all values of s(â) between 0 and 1, inclusively, and implies that the individual spheres of the non-wetting phase remain discrete and non-overlapping. To determine the effective generalized permittivity of a system containing this type of inclusions, the heterogeneous inclusions are replaced by homogeneous pores with the appropriate value of e)(s(1)). Therefore, the dielectric response of such a porous 62 medium is described by (0) * Eg M - 8g 2 Eg+ Cpore (am) g(c0), (e)(s(l)) am)]. (4.20) ml = For the second case, the geometrical configuration of pore fluids within the individual pore spaces is represented by a two component, confocally-layered spheroid. Let the outer spheroidal shell be composed of the wetting phase and the inner spheroid consist of the non-wetting fluid. Stepin (1965) derived an expression for the response of a medium containing identical confocally-layered ellipsoidal inclusions. Expanding this result, the following expression is obtained for a porous medium containing confocally-layered pore spaces: - (am) h(e0),E1)2),sW,am)] + = (4.21) m=l where (i) (2) (1) h e t Â° m 5 g â g â g â 1 (E(Â°) N e(2) 1 ((o) D ), c2) (â) am) + sj), am) Nl(40),e),E2),s),am) m â â (cÂ°), 2 N E(2) (â) (eÂ°), 2 D E(2) (â), am) am) _)) ((2)...e(0)) x = 63 (4 22) (2)_ ( W ) [B(sj),aj 8 Dl(E),e1),e2),sj),am) [A(am) (E(Â°) (2)) [B(sW,am) + 5WA(am) (1 - (4.23) (e2)_e.42) x A(Xmj (e(2) 8(i)) (c(0) c(2)) (e0),E),E2),sl,a)= 8c2)(E)_E0)) 2 N [B(sj), e(2)] sj) - (4.24) 4(1_sj))(c2)_e0))x am) (41)_e2))_(21)+E2))] (4.25) and (4Â°), 2 D e(2) [A(am) (42)_E0)) The term B($j),C) is axis of the - s), am) = [B ((â) (42)+80))] + $i) am) [1 - ((â)_ e(2)) - (c(2)+eO))] (A(Om))2] ((2)_(i)) ((o)_(2)) (4.26) the depolarization coefficient associated with the unique principal inner spheroid. Evaluation of B(s),a J requires the determination of the 11 aspect ratio of the inner spheroid which is a function of both 2). For both types of heterogeneous inclusions, the cases and am (See Appendix =0 and = 1 reduce to the expressions for homogeneous inclusions composed of the wetting phase and non wetting phase, respectively. From the above results, it can be seen that the mathematical formulations for the 64 dielectric and elastic properties of a partially saturated porous rock differ for the various geometrical configurations of the pore fluids within the individual pore spaces. To appreciate the significance of these differences, it is necessary to simulate the physical behavior of a model system using these expressions. In the following analysis, this approach will be used to illustrate the factors that control the relationship between porescale fluid distribution and macroscopic physical properties. Numerical Simulations A partially saturated porous rock containing air and deionized water will be used in the following numerical simulations. In order to employ the preceding mathematical expressions, it is necessary to specify the constituent properties. In the case of electromagnetic wave propagation, it will be assumed that the rock matrix and pore fluids are simple, linear materials where o(j) (j=O,l,2); (4.27) e4i) and ai) are the apparent permittivity and conductivity of the j th component, respectively. Both apparent quantities are taken to be independent of o. This is a reasonable approximation for the dielectric behavior of commonly occurring pore fluids, and both e) and J) are easily obtained. However, due to complex physical and chemical processes taking place at the rock-water interface, determination of the rock matrix parameters requires special consideration. Knight and Endres (1990) showed that 65 the adsorbed water in a water-wetted sandstone should be considered as part of the rock matrix and the parameters used to describe this system should reflect this wetted state. This is an important factor since the parameters describing a wetted and dry rock matrix can differ significantly, particularly in the case of sandstones with large values of surface area to volume ratio (SN). By incorporating the adsorbed water into the rock matrix system, the S used in the theoretical models refers to the level of bulk water saturation (Sbw). The parameters used to describe the dielectric properties of the constituents are given in Table (4.1). The wetted rock matrix values were obtained from an analysis of dielectric data for a sample of Berea sandstone. The frequency used in the modeling computations is 105 kHz. The parameters used for modeling the elastic properties are given in Table (4.2). It should be noted that the difference between dry and wetted matrix elastic parameters was observed by ToksĂ¶z et al. (1976); their value for wetted Berea sandstone is used in this analysis. Initially, the effects of the pore fluid configuration within the individual pore spaces will be investigated. For each type of pore fluid configuration, a series of models containing a single pore geometry (i.e. all inclusions in a given model have the same value of cc) is generated; models within a series are varied by changing the aspect ratio of the pore space. For a given model, the physical property of interest is determined as a function of Sb; for presentation purposes, these results are normalized to the Sb= 0.00 value of the quantity under consideration. It is assumed that the value of local Sb (i.e. the value of Sbw for an individual pore space) for all pore spaces is identical to the overall value of 5 bw when a heterogeneous inclusion type is utilized. The total porosity Ăpore is 66 maintained at 0.01 in all models. Let us first consider the relationship between pore fluid configuration and the dielectric response. Figures (4.1) and (4.2) give the computed values of the apparent dielectric constant (i.e. 1=C/E vacuum) and the apparent conductivitya, respectively, as a function of Sb for models employing homogeneous inclusions. It can be observed that both ic and a display nearly linear variations with changing Sb. Figures (4.3) and (4.4) show the computed values of 1c and a, respectively, as a function of Sbw for models utilizing heterogeneous inclusions filled with an effective fluid consisting of air bubbles in water. These results differ significantly from those employing homogeneous inclusions in that distinctly non-linear relationships are observed for both ic and a. Figures (4.5) and (4.6) illustrate the computed values of 1c and a, respectively, as a function of Sb for models using heterogenous inclusions consisting of an air inner spheroid and confocal water shell. Once again, these results differ from those obtained with homogeneous inclusions in that distinctly non-linear relationships are observed for both ic and a. Further, comparison of results obtained with the two types of heterogeneous inclusions show that the variation in the pore fluid configuration drastically alters the functional dependence on Sb. Another way to alter the pore fluid configuration within the heterogeneous inclusions is to reverse the roles of the wetting and non-wetting phases. The effects of wetting phase variation on the dielectric properties of partially saturated rock has been experimentally observed by Poley et al. (1978). In the above simulations, the water phase was taken to be the wetting phase. While air cannot realistically act as the wetting phase 67 in the pore fluid configurations used in this analysis, its dielectric properties do not significantly differ from those of many hydrocarbons in comparison to the dielectric properties of water. Therefore, it is useful to consider the effect of reversing the role of the two pore fluids used in this study. Figures (4.7), (4.8) and (4.9), (4.10) show the computed values of ic and a, respectively, for the two types of heterogeneous inclusions when air-wetted. Comparison of the air- and water-wetted states for either type of heterogeneous inclusion shows that variation in wetting phase significantly affects the form of the Sbw dependence. It should be noted that in both cases the water-wetted matrix parameters were used; it is expected that rock matrix parameters should vary with changing wetting conditions. However, this effect cannot be quantified at present. Let us now consider the relationship between pore fluid configuration and the elastic behavior of the model. In the derivation of the mathematical formulations, it was shown that the contents of a heterogeneous inclusion can be replaced by an effective pore fluid having elastic moduli that are independent of the geometrical configuration of the pore fluids within the inclusion. Therefore, it is only necessary to compare models containing homogeneous and these âgenericâ heterogeneous inclusions when analyzing the effect of pore fluid configuration on the macroscopic elastic behavior. Figures (4.1 1)-(4. 14) illustrate the computed values of k*, V and V , respectively, as a function of Sb for models containing homogeneous inclusions. Models employing homogeneous inclusions exhibit nearly linear variations of elastic properties with changing Sbw. Figures (4.15)(4.18) show the computed values of k*, V and V , respectively, as a function of SbW for models containing the âgenericâ heterogeneous inclusions. These results differ significantly from those obtained with homogeneous inclusions in that distinctly non 68 linear relationships are obtained between the elastic properties and Sb. In particular, little change occurs in the elastic properties until high levels of Sb (i.e. in excess of 0.95) are obtained; this is when a signficant increase occurs in the bulk modulus of the effective pore fluid. Let us now consider the effects of local variation in saturation conditions on the effective physical properties by using three simple schemes for varying local Sb within the pore geometry spectrum. These simple schemes for increasing overall Sb from 0.00 to 1.00 are as follows: (1) water preferentially fills the largest cc geometries; (2) water preferentially fills the smallest a geometries; and, (3) all geometries fill in unison. In all cases, homogeneous inclusions will be used to represent the individual pore spaces. These schemes are schematically shown in Figure (4.19). The pore spectrum used is an estimate for the pore geometries present within a Berea sandstone based upon Cheng (1978) and is given in Figure (4.20). Figures (4.21)-(4.26) show the computed ic, a, k*, V and V , respectively, as a function of Sb when these three schemes are used to vary the pore-scale fluid distribution. These results illustrate the significant differences which can occur in the functional relation between a given physical property and Sbw when the sequence in which local saturation conditions vary is changed. In particular, changes in the saturation conditions in the smaller aspect ratio geometries have an effect on the dielectric and elastic behavior which is well out of proportion to their percentage of the total porosity (i.e. at high Sb value for scheme 1 and at low Sb values for scheme 2). Unlike the other effective physical properties in this study, both 69 â4 and v are themselves functions of effective quantities. Therefore, variations in both V, and V can be viewed in terms of the relative changes in these quantities. From Equations (4.11) and (4.12), it can be shown that dV dSb = i i 2V = dV dk* dS,w * + dSb + d.L* dV d,.t* dSb a_ di.tâ 3V dSb - dp* + dp* dSbw Ăžrâore V (p(w) - p(a)) 2 (428) and dV dSb = where (w) â dv; dp* dj.t* dSb di.t - Â±.4y dp* dp* dSb ĂžPoreVs (p(w) - p(a))], (4.29) and p (a) are the densities of water and air respectively. In general, , dk* O dSbw (4.30) and (4.31) Therefore, it follows that 70 1 2V â2â di? >0 dSb 3V dSb (4.32) d i 2V: dSb (4.33) â and 0. Thus, the existence of a relative minimum or maximum indicates values of Sb where there is a change in the nature of the dominant mechanism that controls the relationship between V, or v: and Sb. A positive slope indicates that changes in the elastic moduli dominate; a negative slope means that variation in the effective density is controlling parameter. In general, variations in elastic wave velocities due changes in the saturation conditions in the larger cx pores are controlled by changes in paâ. Conversely, changes in the saturation conditions in the smaller cx inclusions cause variations in V and V that are dominated by changes in the effective elastic moduli. Conclusions Inclusion-based formulations describing the elastic and electromagnetic wave propagation in a partially saturated porous rock have been derived by depicting this system as a rock matrix containing inclusions representing the individual pore spaces. By considering various styles of pore-scale fluid distribution, it has been found that drastically different functional relationships between these physical properties and overall saturation level are obtained while maintaining fixed constituent properties and pore 71 structure. Numerical simulations performed with these formulations have illustrated the significance of pore-scale fluid distribution in determining elastic and electromagnetic wave propagation within partially saturated porous rocks. In particular, it has been demonstrated that the pore fluid configuration within the individual pore spaces and the pore geometries in which local saturation conditions are varying are critical factors in determining the physical properties of partially saturated rocks. Further, changes in the saturation conditions in the smaller aspect ratio geometries have an affect on the physical response of the porous medium that is well out of proportion to their percentage of the total porosity. In the next chapter, these basic results will be used to analyze the experimentally measured dielectric and elastic response of partially saturated porous rocks. 72 Rock Matrix Deionized Water Ak 1Ca 7.604 80. 1.01 a(in) 3 1.387x10 7.75x 10 0. Table 4.1: Dielectric constants and conductivities used for the constituents of the model sandstone. k (in -N-) 2 m t (in N) p(in) Table 4.2: Rock Matrix Deionized Water Ak 3.2 x 1010 2.237 X i0 1.550 x 2.45 x 1010 0. 0. 2650. 998. 1.293 Elastic moduli and densities used for the constituents of the model sandstone. 73 1.1 0. . 0.9 N Cu E L. 0 0.7 z 0.5 0.0 0.2 0.4 0.6 0.8 1.0 Bulk Water Saturation Figure 4.1: Computed 1cI1c(Sbw=O.OO) as a function of Sb for a series of simple models containing homogeneous pores. Aspect ratio (cc) of the pore geometry present in a given model is denoted. 74 1.25 Cu E 0) Cl) 001 1.15 0) N 0.10 0.05 1.05 0 z 0.50,1.00 0.95 0.0 âą âą 0.2 âą 0.4 âą 0.6 0.8 1.0 Bulk Water Saturation Figure 4.2: Computed c/(Sb=O.OO) as a function of Sb for a series of simple models containing homogeneous pores. Aspect ratio (a) of the pore geometry present in a given model is denoted. 75 1.3 0. 0. a) N 1.1 0.9 Cu E 0 z 0.7 0.5 0.0 0.2 0.4 0.6 0.8 1.0 Bulk Water Saturation Figure 4.3: Computed 1/ic(Sbw=O.OO) as a function of Sb for a series of simple models containing pores filled with the water-wet effective pore fluid. Aspect ratio (x) of the pore geometry present in a given model is denoted. 76 1.25 Cu E 0) 1.15 a) N Cu E 1.05 0 z 0.95 0.0 0.2 0.4 0.6 0.8 1.0 Bulk Water Saturation Figure 4.4: Computed aIa(Sb=O.OO) as a function of Sbw for a series of simple models containing pores filled with the water-wet effective pore fluid. Aspect ratio (a) of the pore geometry present in a given model is denoted. 77 2.0 0 CD 1.5 N CD E 1.0 L.. 0 z 0.5 0.0 0.2 0.4 0.6 0.8 1.0 Bulk Water Saturation Figure 4.5: Computed 1/1c(Sb=O.OO) as a function of Sb for a series of simple models containing water-wet confocally-layered pores. Aspect ratio (ci) of the pore geometry present in a given model is denoted. 78 1.25 E 0) Cl) 1.15 0 a) N 1.05 0 z 0.95 0.0 0.2 0.4 0.6 0.8 1.0 Bulk Water Saturation Figure 4.6: Computed o/(Sb=O.OO) as a function of Sb for a series of simple models containing water-wet confocally-layered pores. Aspect ratio (ce) of the pore geometry present in a given model is denoted. 79 CD Cu 1.5 N Cu E h. 1.0 0 z 0.5 0.0 0.2 0.4 0.6 0.8 1.0 Bulk Water Saturation Figure 4.7: Computed 1c/1â(Sbw=O.OO) as a function of Sb for a series of simple models containing pores filled with the air-wet effective pore fluid. Aspect ratio (ce) of the pore geometry present in a given model is denoted. 80 Cu 1.5 V a, N CU E 1.0 0 z 0.5 0.0 0.2 0.4 0.6 0.8 1.0 Bulk Water Saturation Figure 4.8: Computed 1/1(Sb=O.OO) as a function of Sb for a series of simple models containing air-wet confocally-layered pores. Aspect ratio (ce) of the pore geometry present in a given model is denoted. 81 Cu E C) Cl) a) N E 0 z 0.0 0.2 0.4 0.6 0.8 1.0 Bulk Water Saturation Figure 4.9: Computed a/a(Sb=O.OO) as a function of Sb for a series of simple models containing pores filled with the air-wet effective pore fluid. Aspect ratio (ct) of the pore geometry present in a given model is denoted. 82 E Cl) N E 0 z 0.95 0.0 0.2 0.4 0.6 0.8 1.0 Bulk Water Saturation Figure 4.10: Computed cIa(Sb=O.OO) as a function of Sbw for a series of simple models containing air-wet confocally-layered pores. Aspect ratio Qx) of the pore geometry present in a given model is denoted. 83 1.9 1.7 0.01 1.5 0 z 0.05 0.10 1.1 0.50,1.00 0.9 0.0 âą 0.2 âą 0.4 âą 0.6 0.8 1.0 Bulk Water Saturation Figure 4.11: Computed k*/k*(Sbw=O.OO) as a function of Sbw for a series of simple models containing homogeneous pores. Aspect ratio (cL) of the pore geometry present in a given model is denoted. 84 1.125 1.100 0.01 1.0751.0500.50 E z 0.10 1.0251.000 - 0.50,1.00 0.9750.0 âą 0.2 âą 0.4 âą 0.6 0.8 1.0 Bulk Water Saturation Figure 4.12: Computed * /*(sbwo 00) as a function of SbW for a series of simple models containing homogeneous pores. Aspect ratio (ct) of the pore geometry present in a given model is denoted. 85 1.20 1.15 > V w N 0.01 1.10 0.05 1.05 0.10 1.00 0.50,1.00 0.95 0.0 âą 0.2 âą 0.4 âą 0.6 0.8 1.0 Bulk Water Saturation Figure 4.13: Computed V/V(Sb=O.OO) as a function of Sb for a series of simple models containing homogeneous pores. Aspect ratio (a) of the pore geometry present in a given model is denoted. 86 1.06 1.05 0.01 1.04 N E 1.02 z 1.01âą 0.10 0.05 1.00 0.50,1.00 0.99 0.0 âą âą 0.2 âą 0.4 âą 0.6 0.8 1.0 Bulk Water Saturation Figure 4.14: Computed v:/v:(sbW=o.oo) as a function of Sb for a series of simple models containing homogeneous pores. Aspect ratio (a) of the pore geometry present in a given model is denoted. 87 1.9 1.7 N 15 0.01 0.9 0.0 0.05,0.10,0.50,1.00 âą 0.2 âą 0.4 â âą 0.6 0.8 1.0 Bulk Water Saturation Figure 4.15: Computed k*/k*(Sbw=O.OO) as a function of Sb for a series of simple models containing heterogeneous pores. Aspect ratio (cL) of the pore geometry present in a given model is denoted. 88 1.125 1.100 1.075 . 0.01 N Co â 1.050 E 0 z 1.025 1.000 ) 0.05,0.10,0.50,1.00 F 0.975 0.0 0.2 0.4 0.6 0.8 1.0 Bulk Water Saturation Figure 4.16: Computed .t*/$t*(Sbw=O.OO) as a function of Sb for a series of simple models containing heterogeneous pores. Aspect ratio (a) of the pore geometry present in a given model is denoted. 89 1.20 1.15 > 0 C) N 1.10 Cu E 1.05 0 z 1.00 0.95 0.0 0.2 0.4 0.6 0.8 1.0 Bulk Water Saturation Figure 4.17: Computed V/V(Sb=O.OO) as a function of Sb for a series of simple models containing heterogeneous pores. Aspect ratio (ce) of the pore geometry present in a given model is denoted. 90 1.06 1.05 Cl) 1.04 V a) N 0.01â> 0.05,0.10,0.50,1.00 0.99 0.0 0.2 âą âą âą âą 0.4 0.6 0.8 1.0 Bulk Water Saturation Figure 4.18: Computed Vâ/V(Sbw=O.OO) as a function of Sbw for a series of simple models containing heterogeneous pores. Aspect ratio () of the pore geometry present in a given model is denoted. 91 Increasing SW E a) U a) a) -cU a) 2 a) U 000 ... 000000 âąâąâąâąâąâą âąâąâąâąâąâą 000000000 âąâąâą 000000 âąâąâąâąâąâą 000 âąoo âąâąoâąâąâą 000 âą ooâą âą oâąâąâą c cz cz czD czD cZD c ED SW Figure 4.19: = cz z cz - cz z 0.0 cz cz cz - - â - cz SW= 1.0 Schematic diagram of the three fluid distribution schemes using homogeneous pores for varying local Sb level within the pore geometry spectrum. Pore space is illustrated with three different cc spheroids. Shading=water white=air. Scheme one-water preferentially fills large cc pores; scheme two- water preferentially fills small cc pores; scheme three-water fills all pore geometries in unison. 92 1W 1 d 10 Oc =â:- cu a E 1O-ââą . . . iĂŒ . . o_ io âą . - 6. 10:5 03W2 Aspect 0Â° Ratio Figure 4.20: Pore geometry spectrum in volume fraction for model sandstone. 93 30 25 Cu 0. 0. 20 Cu 15 10 5 0.0 0.2 0.4 0.6 0.8 1.0 Bulk Water Saturation Figure 4.21: Computed w as a function of Sbw for the model sandstone using the three fluid distribution schemes for varying local Sbw level. 94 1 .30e-3 1 .20e-3 c 1.lOe-3 E 0) Cl) 1 .OOe-3 9.OOe-4 0.0 0.2 0.4 0.6 0.8 1.0 Sw Figure 4.22: Computed c as a function of Sb for the model sandstone using the three fluid distribution schemes for varying local Sb level. 95 2.20+10 2.OOe+1 0 1.80e+10 E 1.60e+10 . 1.40eĂ·10 1.20e+10 1.OOe.i-10 0.0 0.2 0.4 0.6 0.8 1.0 Bulk Water Saturation Figure 4.23: Computed k* as a function of Sb for the model sandstone using the three fluid distribution schemes for varying local Sbw level. 96 1.40e+10 / C1 E Scheme #2 Scheme #3 1.30e+10 z 1.20e+1O Scheme #1 1.lOe+10 0.0 âą â 0.2 âą âą 0.4 â 0.6 âą 0.8 1.0 Bulk Water Saturation Figure 4.24: Computed tâ as a function of SbW for the model sandstone using the three fluid distribution schemes for varying local Sb level. 97 4200 4000 Cl) 3800 0. >1 3600 3400 0.0 0.2 0.4 0.6 0.8 1.0 Bulk Water Saturation Figure 4.25: Computed V as a function of SbW for the model sandstone using the three fluid distribution schemes for varying local SbW leveL 98 2500 Scheme #2 2400- E Cl) > Scheme 2300- Scheme #1 22000.0 âą âą 0.2 âą 0.4 âą 0.6 0.8 1.0 Bulk Water Saturation Figure 4.26: Computed V as a function of Sb for the model sandstone using the three fluid distribution schemes for varying local Sbw level. 99 Chapter Five The Effects of Pore-Scale Fluid Distribution on Elastic and Electromagnetic Wave Propagation in Partially Saturated Porous Rocks: An Analysis of Experimental Results Introduction A number of experimental studies have clearly demonstrated the dependence of macroscopic physical properties on pore-scale fluid distribution (Domenico, 1976, 1977; Bourbie and Zinszner, 1984; Knight and Nur, 1987a; Longeron et al., 1989; Knight and Nolen-Hoeksema, 1990); differing styles of pore-scale fluid distribution were obtained in these laboratory experiments by changing the saturation technique. The theoretical importance of pore-scale fluid distribution in determining the dielectric and elastic response of a partially saturated porous rocks was demonstrated in the last chapter. The mathematical formulation employed in that analysis will now be used to study the experimentally measured elastic wave velocities and apparent dielectric constant of partially saturated porous materials. Of particular interest are the results of Knight and Nur (1987a) and Knight and Nolen Hoeksema (1990) that give the measured apparent dielectric constant and elastic wave velocities, respectively, for similar samples of a tight gas sandstone undergoing the same two saturation techniques. Before the existence of these results, individual experimental studies have been concerned with the effect of differing saturation techniques on a single 100 macroscopic physical property. Further, a variety of different saturation procedures and porous media samples were used in these studies; hence, there is difficulty in comparing the results of these experiments. The data set obtained by Knight and co-workers permits a direct assessment of the variation in the dielectric and elastic response of the same rock type undergoing the same saturation methods. Further analysis of the relationship between elastic wave velocities and differing styles of pore-scale fluid distribution is undertaken by considering the experimental data obtained by Domenico (1977). In his experiment, Domenico subjected a glass bead packing to two saturation techniques that differ from those employed in the tight gas sandstone experiments. It is shown in this study that the physical properties of both the tight gas sandstone and the glass bead packing are replicated by utilizing simple models to describe the geometrical configuration of the constituents. In addition, it is seen that a single geometrical model for pore-scale fluid distribution can be used simultaneously to estimate the dielectric and elastic behavior of the tight gas sandstone. Tight Gas Sandstone: Experimental Results Knight and Nur (1987a) measured the apparent dielectric constant (iĂ§) of a tight gas sandstone sample from the Spirit River Formation in the Alberta Basin; Knight and Nolen Hoeksema (1990) measured the compressional wave velocity (Vp) and shear wave velocity (Vs) on a similar sample from the same formation. The experimental procedures employed are described in detail in these two studies. The sample used in the dielectric measurements has a porosity of 0.070 and a gas permeability of 7.18 jt D; the sample used in the study of the elastic wave velocities has a porosity of 0.052 and a permeability of 1 101 i D. A complete petrographic description of the Spirit River samples is found in Walls (1982). In both laboratory studies, measurements were made during a continuous imbibition drainage experiment (where imbibition is the displacement of air by water in a water-wet system and drainage is the displacement of water by air). Starting with a desiccator dried sample, S was initially increased through imbibition by adsorption of water vapor, followed by the soaking of the sample in deionized water. The maximum value of S obtained through this procedure was approximately S 0.90 in both experiments. Water saturation was subsequently decreased through drainage by evaporative drying. As S was varied, its value was determined by means of sample weight. Figures (5. 1)-(5.3) display the measurements of ica, VP and V, respectively, obtained during the imbibition-drainage cycle. The dielectric data were measured at a frequency of 105 kHz; the central frequencies for the V and V measurements were 1 MFIz and 600 kHz, respectively. Each data set clearly illustrates the dependence of the functional relationship between a physical property and SW on the saturation technique used. There are clearly two distinct regions in both the apparent dielectric constant and elastic wave velocity data sets that can be related to the location and nature of the water present in the pore space. Knight and Nur (1987a,b) defmed a critical value of SW (SWO) separating these two regions in the dielectric response (i.e. SÂ° = 0.15); these two regions are labelled in Figure (5.1). In region 1 (S <SWÂ°), the water present in the rock exists as monolayers of adsorbed water coating the surfaces of the pore space. In region 2 (SW > SW ), 0 a bulk water phase filling the central portion of the pore space is also present. For the sample used in the dielectric measurements, SWÂ° has been shown to correspond to a surface layer 102 of bound water, approximately 1 nm in thickness (Knight and Nur, 1987a,b). The same distinction can be made in an analysis of the elastic data shown in Figures (5.2) and (5.3). In this case, the value of SÂ° cannot be precisely defmed and a zone of transition in which S occurs is denoted (i.e. 0.25 Sj 0.35). In region 1, where only bound water is present, a significant increase in iCa and a decrease in the elastic wave velocities, particularly V (Clark et al., 1980), occurs as S increases towards SÂ°. Variations in the physical properties within this saturation range are due to complex rock-bound water interactions. Since rocks in situ usually contain the bound water phase, the physical behavior of rocks in this region of S values will not be considered in this study. Of interest in this study is the theoretical analysis of the hysteresis that occurs in both data sets at higher S values. This hysteresis in the physical properties has been attributed to changes in the geometrical distribution of fluids within the pore space that occur during the imbibition-drainage cycle (Knight and Nur, 1987b; Knight and Nolen-Hoeksema, 1990). Based on earlier works by Haines (1930) and Foster (1932), the following model was proposed by Knight and Nur (1987a) in order to explain the results of their dielectric measurements. At low levels of S, as water content is increased by the adsorption of water vapor, the water preferentially adsorbs on the surfaces of the pore space. This adsorption leads to the development of thick surface layers of water along the outer portion of the pore volume with a thin, interconnected gas phase within the central portion. This arrangement of gas and water means that the individual pore spaces are in a partially saturated state. When the imbibition process is continued at higher levels of SW by soaking the sample, the bulk water phase tends to maintain this established geometrical configuration by continuing to move 103 along and coat the surfaces. However, this geometrical configuration is a metastable state persisting until closure occurs between the surface layers of bulk water. At this point, the rearrangement of the bulk water and air to a more thermodynamically stable configuration occurs, causing the break-up of the interconnected gas phase and the filling of the central pore volumes with the bulk water phase. The model proposed for the drainage process is quite different. During drainage, the bulk water empties from each individual pore space as it is accessed, removing all but the first few tightly-held monolayers of bound water. Hence, the geometrical configuration of the pore fluids that occurs during imbibition is not re-created during drainage. The imbibition process favors the development of pores with co-existing gas and bulk water phases; the drainage process tends to develop either fully bulk water-saturated or fully gassaturated pore spaces. Since the object of this study is to assess the effects associated with the geometrical distribution of the bulk water phase, we are only concerned with the behavior of the tight gas sandstone in the range of the S values between 50 and 1.00. In this region, the discussion will be in terms of bulk water saturation (Sb) over a range from 0.00 to 1.00. Further, the bound water phase is considered to be part of a wetted rock matrix system (Knight and Endres, 1990); this implies that the total porosity of the samples must also be corrected for the existence of the bound water phase. If the value of SWÂ° is taken to be 0.15 and 0.35 for the dielectric and the elastic measurements, respectively, then the corrected porosity (Ă, the total porosity occupied by the bulk water phase) for the dielectric and elastic samples is 0.0595 and 0.0338, respectively. 104 Tight Gas Sandstone: Numerical Simulation Let us first consider the drainage process. Given that the segregation of pore fluids occurs, homogeneous inclusions composed of either air or bulk water are used to model the pore-scale fluid distribution. The effective permittivity and elastic moduli of such a system with homogeneous pore spaces is calculated by using Equation (4.4) and Equations (4.7) and (4.8), respectively. The overall value of Sb is varied by simultaneously changing the local Sb level in all pore geometries such that emptying occurs in unison across the entire pore geometry spectrum. This scenario for the pore-scale fluid distribution is schematically shown in Figure (5.4). Given that the topology of the pore space within a rock necessitates essentially simultaneous drainage of all aspect ratio geometries in order to access the total pore space, this scheme is considered to be realistic for modeling the experimental data. For the imbibition process, it is necessary to treat each individual pore as being simultaneously occupied by air and bulk water until the rearrangement of pore fluids occurs. Given the configuration of thick surface layers of bulk water and an interconnected air phase in the central volume of the pores, confocally-layered spheroids consisting of a bulk water outer shell and an inner air spheroid are reasonable inclusion geometries to use in simulating the process. Further, it will be assumed that all pores with a given shape will have identical levels of local Sb prior to the critical rearrangement phenomenon. The effective permittivity for a system containing confocally-layered pore spaces is computed from Equation (4.21). Since the pore spaces are heterogeneous, Equations (4.16) through (4.18) are used to determine the effective elastic moduli. To simulate the pore-scale fluid distribution which occurs during the imbibition process prior to the critical rearrangement of pore fluids, the following idealized geometrical 105 arrangement of pore fluids was employed; this geometrical configuration is schematically illustrated in Figure (5.5). As Sbw starts to increase from 0.00, the lower aspect ratio geometries, which are generally associated with the pore throats (in the case of the tight sandstone model, this will be assumed to be a 0.5), are initially filled to some level of local Sbw. The local levels of Sbw specified for the various pore throat shapes in this simulation are given in Table (5.1). This establishes the presence of the narrow gas pockets within the pore throats. This geometrical configuration within the pore throats is maintained as the overall level of Sbw of the medium is increased by confocally filling the remaining pore geometries until S the point at which the rearrangement of pore fluids , occurs, is achieved. The critical value of overall water saturation S is taken to be 0.65 and 0.80 in the case of the dielectric and velocity models, respectively. After the rearrangement of the pore fluids, it will be assumed that the pore-scale fluid distribution is identical to that used to model the drainage process at the given value of Sb. The constituent properties used for modeling the dielectric and elastic properties of the tight gas sandstone are given in Tables (5.2) and (5.3), respectively. The frequency used in the computation of the dielectric properties is 105 kHz. Figure (5.6) illustrates the pore geometry spectra used to describe the two samples. These spectra are based on the work of ToksĂ¶z et aL (1976) and Cheng and ToksĂ¶z (1979). Since different samples were used in each of the experiments, small differences in the fractional volume for a given aspect ratio are permitted between the pore geometry spectra representing the two samples. In both cases, the total porosity used is corrected for the bound water being incorporated as part of the rock matrix system. 106 Figure (5.7) shows the computed apparent dielectric constant (ic); Figures (5.8) and (5.9) show the computed compressional and shear wave velocities (â4 and V), respectively. In all cases, the experimental data is superimposed on the model response. While there are differences between the theoretical models and experimental data, the basic form of both are in good agreement. The model accurately predicts the functional form and magnitude of the hysteresis which occurs in both the apparent dielectric constant and elastic wave velocities as the sample undergoes the imbibition drainage cycle. The differences - which are observed can be attributed to inaccuracies both in the geometrical schemes used to represent the pore-scale fluid distribution and in the parameters used for the constituents and pore geometry spectra. However, it is significant that the relatively simple scenarios used to depict the geometrical configurations occurring within the rock give good results. In particular, it can be observed that the existence of the thin interconnected air phase throughout the pore space during the imbibition technique prior to the rearrangement of the pore fluids causes a significant enhancement in Ka and a depressed, almost constant value of V, and V. Further, the drainage technique gives a much more uniform variation in physical properties with changes in SbW. Glass Bead Packing: Experimental Data Domenico (1976, 1977) investigated the functional relationship between elastic wave velocities and S for a glass bead packing by employing two different saturation methods that he referred to as the flow and imbibition techniques. The experimental procedure used is given in detail in the above references and is briefly summarized below. In the flow technique, the sample was initially saturated with nitrogen gas; the gas was subsequently 107 totally displaced by a brine solution. Partial saturation was achieved by injecting a series of gas-brine mixtures into the sample. The gas and brine were mixed prior to injection; however, the two pore fluids remained immiscible during the emplacement procedure. Approximately ten pore volumes of the prepared fluid mixture were passed through the glass bead packing between measurements. The imbibition method initially utilized the flow method to obtain a partially saturated state. Subsequently, the fluid exit valve was closed and pure brine was injected into the glass bead packing until the gas phase present within the pore space was completely dissolved into solution. This results in a significantly increased pore fluid pressure (pf) within the porous medium. To obtain a partially saturated medium, p was decreased and gas was exsolved from solution. During this process, the confining pressure (Pc) was adjusted such that the effective pressure (Pe=Pc-Pf) was constant for all measurements. During both saturation techniques, the level and uniformity of S within the glass bead packing was monitored by means of an X-ray system. Figures (5.10) (5.13) illustrate the experimental results for the elastic wave velocities - as a function of S; these measurements were performed at a central frequency of 200 kHz. Over the range of S values between 0.70 and 1.00, the two saturation techniques gave distinctly different results for V at both effective pressures. Although the flow technique produced a less uniform distribution of S values along the length of the sample, these variations are not sufficient to explain the observed differences between the two saturation methods. For S values between 0.88 and 1.00, the flow technique gives a very uniform S distribution while differences between the two methods are the most 108 pronounced. There are indications in the V measurements that the two saturation techniques are different; however, the data are not conclusive. It is hypothesized in this study that Domenicoâs flow technique results in a pore-scale fluid distribution similar to that obtained during the imbibition method employed by Knight and co-workers. By initially saturating the medium completely with brine, it is established as the wetting phase throughout the glass bead packing. Subsequent injection of a distinct gas phase results in it being placed in the central portion of the pore spaces. Up to some critical level of S, (Se), the pressure within the gas can overcome other thermodynamic factors in order to establish and maintain an interconnected gas phase. Once again, such a configuration of fluids is thermodynamically metastable; and, gas pressure cannot maintain gas interconnectivity beyond S. However, there is a significant difference between these two saturation techniques. Knightâs method results in a continuous sequence of changes in the pore-scale distribution of fluids which represents a single realization of the metastable process. Hence, physical properties display a clear functional relationship as S is varied. However, this metastable process is sensitive to the saturation history of the porous medium. By flushing large volumes of fluids through the porous medium, Domenicoâs technique produces relatively independent realizations of the metastable process for each measurement. Hence, this technique can produce significantly different pore-scale fluid distributions for similar values of S. In particular, the value of S at which S occurs varies between realizations of this process. This is reflected in the scatter observed in the elastic wave velocity data and the variation observed in local S (see Figure 22 in Domenico,1977) over the range of S values between 0.70 and 0.90. 109 Domenicoâs imbibition technique results in a style of pore-scale fluid distribution which is significantly different from those previous discussed. The exsolution of a gas phase from solution to obtain high values of S should result in the formation of individual gas bubbles throughout the pore space. Further, the creation of the gas bubbles should be relatively uniform through the pore structure of the rock. Glass Bead Packing: Numerical Simulation The pore-scale fluid distribution resulting from Domenicoâs flow method is simulated by means of the same geometrical scheme used to represent Knightâs imbibition technique. The generation of gas bubbles implies that Domenicoâs imbibition technique is characterized by heterogeneous pore spaces. Hence, Equations (4.16) through (4.18) describe the effective elastic moduli resulting from this saturation technique. Further, it will be assumed that the exsolution process produces a uniform distribution of gas bubbles such that the local S within a given pore space is the same as the overall value of S. The elastic properties of the constituents are given in Table (5.4); Figure (5.14) gives the pore geometry spectrum used to represent the glass bead packing. It is assumed that the relative amount of the bound water phase is small in comparison to that of the bulk water phase within the glass bead packing. Therefore, corrections for the volume occupied by the bound water phase were not made; and, the results of the numerical simulation are expressed in terms of total water saturation. In simulating the flow technique, S was taken to be 0.80; and, pore throat geometries were assumed to have aspect ratios of 0.5 or less. 110 Figure (5.15) and (5.16) show the computed compressional and shear wave velocities (V and V), respectively, for the glass bead packing model. The functional relationships between VP and S obtained from the numerical simulation for the two saturation techniques is in good qualitative agreement with the experimental data. In particular, the dominance of the density effect over the lower S, range during the flow method and the dramatic decrease in V as S decreases from full water saturation during the imbibition technique are accurately reproduced. Numerical results for V show that the expected difference between the two saturation methods is very small. Hence, it is not surprising that the experimental data is not conclusive. Conclusions The processes involved in the filling and emptying of a porous rock with fluids are undoubtedly complex. In this study, simple models are employed to describe the porescale fluid distribution resulting from the use of a given saturation technique. The treatment of these saturation methods has been simplified to consider only the basic geometrical elements characterizing the pore space and the contained fluids. The good agreement that is obtained between the theoretical models and the experimental data implies that these geometrical effects are critical factors in determining the physical properties of partially saturated rocks. Further, it is significant that a single geometrical model for the pore-scale fluid distribution was used simultaneously to estimate the elastic wave velocities and apparent dielectric constant of a partially saturated rock. In the theoretical formulations used, it was assumed that communication between individual pore spaces is minimal. Hence, intrapore phenomena control both the elastic and 111 electromagnetic wave propagation in the porous medium. This is achieved experimentally by taldng measurements at sufficiently high sampling frequencies and/or using porous media that possess adequately low permeability. In the cases of the tight gas sandstone and the glass bead packing, these requirements appear to have been satisfied. However, a review of the experimental work in the literature on elastic wave velocities in partially saturated rocks suggests that these conditions commonly do not hold. Drainage by evaporative drying has been used to control overall S in a number of studies (Wylie et al., 1956; Gregory, 1976; Murphy, 1982b, 1984, 1985). The theoretical results for the tight gas sandstone (Figures 5.8 and 5.9) and Berea sandstone (Figure 4.25 and 4.26) imply that elastic wave velocities should vary in a nearly linear manner with changing S. Results obtain by Murphy (1982b, 1984, 1985) show that significant departure from this quasi-linear relationship occurs as sampling frequency decreases. In the case of Massilon sandstone (Murphy, 1982b), the combination of very high permeability (737 mD) and sampling frequencies in the vicinity of 500 Hz gives results that are well described by the phenomenological formulations of Gassmann (195 ib) and Biot (1956) where complete pressure communication within the pore fluid is assumed. However, the behavior of tight gas sandstones (Murphy, 1984) and granite (Murphy, 1985) measured at sampling frequencies of approximately 5 kHz display an intermediate type of functional relationship that cannot be adequately described by the inclusion-based formulations used in this study or the phenomenological approach. Further, there appear to be no adequate methods to determine a priori whether a particular formulation is valid for a given porous rock for a specific set of conditions. 112 Table 5.1: Range of Local Level of S bw cc>10 1 5x10 0.98 2 10â>cc10 0.75 a 10> 2 0.20 4 c 10 1O> 3 0.10 Level of local Sbw maintained in the pore throat geometries during the imbibition process prior to the rearrangement phenomenon. Rock Matrix Deionized Water Ak Ka 13.75 80. 1.01 Ga(inj) 3 2.53x10 5 5.x10 0. Table 5.2: Dielectric constants and conductivities used for the constituents of the tight gas sandstone model. 113 Rock Matrix Deionized Water Air k (in-N-) 2 m 0 2.9x10â 9 2.237x10 5 1.550x10 t(in-N-) 0 2.6x10â 0. 0. p (in5) 2650. 998. 1.293 Table 5.3: Elastic moduli and densities used for the constituents of the tight gas sandstone model. k (in Nâ1 I 2 m t (in -N..) p(inE-) Table 5.4: Glass Bead Packing Brine 2.137 x 1010 2.098 X i0 1.378 x 106 3.053 x iO 0. 0. 2420. 1080. 11.67 Elastic moduli and densities used for the constituents of the glass bead packing model. 114 45. oo0o I (1) 35 : (2) 0 0 â0 A a) A A Â°A0A30 A 25 0. 0 CU 15 L 1 0 : 0 0 0 A I 5 0.0 âą 0.2 Total Figure 5.1: Measured Imbibition Drainage âą 0.4 Water âą 0.6 0.8 1.0 Saturation dielectric COflStaflt (1(a) of a tight gas sandstone undergoing imbibition-drainage cycle (sampling frequency=lO5kHz) (from Knight and Nur, 1987a). Region (1): only bound water phase present; region (2): bulk water present. apparent 115 4500 0 A I I I 4000 I I I A (1) hâ, Cl) 1mb jb ition Drainage (2) AAO H 0 I I I > I I I 3000 D 25000.0 0 âą 0.2 Total 0.4 Water 0.6 0.8 1.0 Saturation Figure 5.2: Measured compressional wave velocity (Vp) of a tight gas sandstone undergoing imbibition-drainage cycle (from Knight and Nolen-Hoeksema, 1990). Region (1): only bound water phase present; region (2): bulk water present. 116 2500 â â I I I I âI I I I 2400 AA : : (â AA (2) :â* Cl) A) - E Imbibition Drainage I I I I II (1) 0 A nâioi.iâQ o0 0 ) Cl) I i I I i 1 I I I AA A >1 2200 2100 0.0 AA â A A %A9po ca 1 000 0 0 :A1A oÂ°Cb I I âą 0.2 âą 0.4 âą 0.6 0.8 1.0 Total Water Saturation Figure 5.3: Measured shear wave velocity (Vs) of a tight gas sandstone undergoing imbibition-drainage cycle (from Knight and Nolen-Hoeksema, 1990). Region (1): only bound water phase present; region (2): bulk water present. 117 00 II . CD I, CD I-. ciq â-Ă·0 c,,CD CD c, C) II CD -io CD . crq C) . 0 0 II 0 U) 0 0 II Cl) azi r:r:csi cicc::i ooc ooc ooc rrrzi ic.::::::.jc::::::::: I (b) â 1 o (c) Figure 5.5: Schematic diagram of the scenario used to simulate the pore-scale fluid distribution during the imbibition process: (a) initial filling of the pore throats; (b) filling of remaining pore geomethes to the critical overall SbW value; (c) all pore geometries filling in unison after rearrangement phenomena. Black=water, white=air. 119 2. 0 A Dielectric Model Elastic Model C AAAA aoooooooci Co OcU -â .1-I CD I0 A A0A o A Câ y 00 -I 1 0 - Aspect 2 - 1 00 Ratio Figure 5.6: Pore geometry spectra (in volume fraction) used for dielectric and elastic models of the tight gas sandstone. 120 50 40 0. 30 20 0.0 0.2 0.4 0.6 0.8 1.0 Bulk Water Saturation Figure 5.7: Computed apparent dielectric constant as a function of bulk water saturation (Sb) for the tight gas sandstone model using the imbibition drainage scenarios. Experimental data are superimposed on the theoretical results. 121 (ic) 4500 4000 Co >1 3000 2500 0.0 0.2 0.4 0.6 0.8 1.0 Bulk Water Saturation Figure 5.8: Computed compressional wave velocity (V) as a function of bulk water saturation (Sb) for the tight gas sandstone model using the imbibition drainage scenarios. Experimental data are superimposed on the theoretical results. 122 2600 2500 Cl) 2400 U) 2300 2200 2100 0.0 0.2 0.4 0.6 0.8 1.0 Bulk Water Saturation Figure 5.9: Computed shear wave velocity (V ) as a function of bulk water saturation (Sbw) for the tight gas sandstone model using the imbibition-drainage scenarios. Experimental data are superimposed on the theoretical results. 123 2200 0 0 Flow A Imbibition 0 00 1900 C,) 1600 ocb 0 1300 Q0 1000W 0.00 0.25 0.50 0.75 A 1.00 Total Water Saturation Figure 5.10: Measured compressional wave velocity (Vp) of a glass bead packing undergoing flow and imbibition methods (PeâÂ°. 33 MPa) (from Domenico, 1977). 124 900 o A Flow Imbibition 800 C,) 9 6 Cl) 0 700 00 0 A 600 0.00 0.25 0.50 0.75 1.00 Total Water Saturation Figure 5.11: Measured shear wave velocity (Vs) of a glass bead packing undergoing flow and imbibition methods (pe=lO. 33 MPa) (from Domenico, 1977). 125 2300 0 A Flow Imbibition 0 2000 Cl) > 1700 QA 0 00 oc 1400 0.00 8* A 0 0.25 0.50 0 0.75 1.00 Total Water Saturation Figure 5.12: Measured compressional wave velocity (Vu) of a glass bead packing undergoing flow and imbibition methods (Pe= 31 .00 MPa) (from Domenico, 1977). 126 1200 o A Flow Imbibition 1100 U) 0 (I) > A 1000 A AA 0 0Â° 0 cP 0 0 900 0.00 0.25 0.50 0.75 1.00 Total Water Saturation Figure 5.13: Measured shear wave velocity (Vs) of a glass bead packing undergoing flow and imbibition methods (Pe 31 .00 MPa) (from Domenico, 1977). 127 100 - a) wE oz 101 () . 3 io- 00 o 10âI ......... âąI ........I 102 Aspect 101 100 Ratio Figure 5.14: Pore geometry spectrum (in volume fraction) used for elastic models of the glass bead packing. 128 2200 2000 (I, 1800 Imbibition 1600 Flow 1400 0.0 âą 0.2 âą 0.4 âą 0.6 0.8 1.0 Total Water Saturation Figure 5.15: Computed compressional wave velocity (V,) as a function of total water saturation (Sw) for the glass bead packing model using scenarios for Domenicoâs imbibition and flow methods. 129 900 880 U) 840 > 820 800 0.0 0.2 0.4 0.6 0.8 1.0 Total Water Saturation Figure 5.16: Computed shear wave velocity (V ) as a function of total water saturation (Sw) for the glass bead packing model using scenarios for Domenicoâs imbibition and flow methods. 130 Chapter Six A Model for Incorporating the Effects of Surface Phenomena in the Dielectric Response of Heterogeneous Media Introduction The presence of surface phenomena can significantly affect the dielectric response of a heterogeneous medium. One method for incorporating the effect of surface phenomena in a theoretical formulation is to redefine the bulk components. In particular, the adsorbed water phase present within a water-wetted sandstone was considered to be part of the rock matrix system; and, the parameters used to describe this system reflect this wetted state (Knight and Endres, 1990). However, this approach does not preserve the geometrical structure of the system; the effects of surface phenomena are uniformly distributed throughout the rock matrix. In order to retain the structure of the system, it is necessary to employ a mathematical formulation that localizes these mechanisms at the interfaces within the system. The region in which surface phenomena occur is generally very thin in comparison to the dimensions of the individual pore spaces. Thus, the effects of surface phenomena can be incorporated into the analysis of an individual inclusion by specifying the appropriate boundary conditions. This approach was first used by OâKonski (1960) in order to introduce surface currents and an associated surface conduction into an analysis of a spherical colloidal particle. Subsequently, a number of studies (Schwarz, 1962; Schurr, 131 1964; Dukhin, 1971; Dukhin and Shilov, 1974; Fixman, 1980, 1983) have employed a variety of different boundary conditions in order to incorporate surface phenomena into the analysis of a spherical particle. To adequately describe the geometry present within a porous medium, it is necessary to consider more versatile inclusion shapes, such as ellipsoids. However, analytical solutions for ellipsoidal inclusion geometries employing the necessary boundary conditions have not been found. OâKonski (1960) offered an approximate solution for the ellipsoidal case based on his results for a spherical particle. He assumed that the effects of the surface phenomena could by represented by an equivalent volumetric mechanism which is uniformly distributed throughout the ellipsoid. Again, this approach radically alters the geometrical configuration of the system under consideration. In this chapter, a new approximation for the case of an ellipsoidal inclusion influenced by surface phenomena is present. Motivated by the technique used by Miles and Robertson (1932) in an analysis of the spherical case, the starting point of this model is a confocally layered effipsoidal system. The outer shell in which the surface phenomena take place is taken to be infinitesimally thin while the number of unit charge carriers and unit dipoles within this region is conserved. In the limit, this results in a homogeneous (i.e. single component) inclusion possessing a surface permittivity representing the surface phenomena. For spherical geometries, this method gives a result identical to that obtained rigorously by means of boundary conditions. Response of a Homogeneous Inclusion with Surface Phenomena Determination of the macroscopic dielectric properties of a heterogeneous medium by 132 means of an inclusion-based formulation requires an analysis of the response of the individual inclusions. The analysis undertaken to this point has employed boundary conditions that exclude surface phenomena at the interfaces between the inclusions and background material. It is possible to incorporate the effects of surface phenomena into this analysis by specifying the appropriate boundary conditions; however, an analytical solution has not been determined for the case of a general effipsoidal inclusion. An approximation for a homogeneous ellipsoidal inclusion when surface phenomena occur can be obtain by initially representing this system as a two-component confocally layered ellipsoidal inclusion. The outer and inner ellipsoids of this confocal system have principal semi-axis lengths of (aĂ§â ),aĂ§â)) and (aĂ§2),$2),aĂ§2)), respectively. The relationship between the semi-axis lengths of the two ellipsoids is given by (a1))2= (a2))2+ where (6.1) (j=l, 2, 3) is the ellipsoidal coordinate associated with confocal ellipsoids. The shape of the outer and inner ellipsoids are uniquely determined by the aspect ratio pairs 2),c42)), 1 (a( (c41),a) and respectively. For a heterogeneous system containing this type of inhomogeneity, Equation (3.6) can be rewritten as VEg+ J (Eg(X) eg(Â°))EO)(x)dv Vmc where the perturbation term now involves the volume integration over the confocally 133 (6.2) layered inclusions. Since the outer shell and inner ellipsoid are composed of homogeneous materials with generalized permittivities 4) and 42, respectively, the evaluation of the perturbation term only requires the determination of the electric field strength within each component. As was done in Chapter Three, E (1)(X) within an individual inclusion is estimated by treating it as an isolated inclusion embedded in an infinite homogeneous background with generalized permittivity subjected to a uniform applied electric field with strength Ea. Stepin (1965) determined the electric field strength within such a system; his analysis showed that the electric field strength within each layer of the inclusion can be decomposed into âhomogeneousâ and âheterogeneousâ components. Further, it was demonstrated that integration of the contribution due to the âheterogeneousâ component over the volume of a given layer is zero. Hence, only the âhomogeneousâ component need be considered in the evaluation of the perturbation term. The âhomogeneousâ component of the electric field strength is 3 0 Ea = - (E(1)e(2)a(2) â cos O 2 ),$ ),4 ),a(1),a(2)) (6.3) eI (6.4) â Ea J=i TIO(4 â within the outer ellipsoidal shell and 0 Ej = Ea E b j=1 [io(4 ),4 ),4 1 134 2 ),a(1),a(2)) C05 within the inner ellipsoid. The terms qo and Thj are defined as follows: loj(ee(1)e(2)a(1)a(2)) = {[A2)(czc2)c42)) _i] !j A2)(xc2)$2))} - {[A1(cc1)a1)) -11 :: A1)(c1$1))} + [1_A1)(xĂ§1),c41))] (â_ )i)i) Aâ)(aĂ§â),aâ)) x )(â- (6.5) and nij(41),42),a(2)) = [A2) ((2)$2)) A ((2)c42)) _1] , (6.6) - where Aâ (c4â ),c4â)) and A2)(aĂ§2),42)) are the depolarization coefficients with respect to the j th axis for the outer and inner ellipsoid, respectively. Using Equations (6.3) and (6.4), it can be shown that the expected value of the volume integral for a confocally-layered inclusion that is oriented in a unifomly random manner is given by KI (cg(x) e) EOkx)dv) - vinc 135 Eae(ac1)41)ac1)) I Nj),$1),$2),a(1),a(2)fl , = Rjfr(),41),42),a(1),a(2)) j. (6.7) J where 1) V(2)(a(2)) N(4Â°)âgâgâ Eâ (2) a(1),a(2)) = (e(2) (O)) {[A2) (a(2) c$2)) - i) 1] - A2) (2)a2) 42) - Vâ - jO)) x v(â) (a(1),a2)} (6.8) and R(e(b) (i) (2) a(1),a(2)) gâgâgâ 1] ) â = {[A2) ((2)2))... i) A2) ((2)c$2)) E(2)} x i] A1)(aĂ§1x1))41)}(aĂ§1) Ă§o aĂ§â)) + A1)(1),x1)) < â [1 A1)( - 1x1))] (e$â) The volume of the outer confocal shell y(1) (a(1),a(2)) = - )(aĂ§ 1 ) 2 4 4â))(e(gâ) ) â ) 1 aĂ§ (6.9) and inner effipsoid are it [(aĂ§ 40 aSâ)) and 136 - ) a5) a52))] 2 (aĂ§ (6.10) V(2) (a(2)) (6.11) (a) aj aĂ§ )), 2 it = respectively. To obtain an estimate for a homogeneous inclusion possessing surface phenomena, it is necessary to consider the limiting case as â* 0. Using a Taylor series expansion, the following approximations for quantities pertaining to the outer ellipsoidal shell can be obtained when is very small relative to a2) (j = 1, 2, 3): (j = )a)aĂ§ aĂ§ ) 1 =2 )a)aĂ§ aĂ§ ) + v(â) (a(1),a(2)) - 1,2,3), C(a(2)) (6.12) (6.13) , 2 (a(2))] it[C (a(2)) + D , (6.14) and Aâ)(Ă§â),c$â))= A2)(cz2),c$2)) + F (a(2)) (j = 1,2,3), (6.15) where C(a(2)) i ( a(2) a(2) a(2) a(2) 223 + 137 ) 2 aĂ§ a(2) a(2) + a(2) )â (6.16) a(2) (2) a 2 + a2)a + aĂ§ aĂ§ 2 (2) a 3 j aĂ§ aĂ§ 2 2 (2) (2) a a 2 3 (2) (2) (2) a a3 + a(2) 1 2 1 a )) 2 (ac 3 (a2))3 (a)) ,( (617) and {L Fi(a(2)) {3{s + (a] + { + (a] + (ac2] 5/2 { (a22] 1 ) 42) $2)) 2 (ac 2 A , ) ac (ac - [s + (aĂ§2] + 3/2 (a] [ + (aS2))] + [ + ($2] + ds 3 (ac2)y] [s + } (4]} (6.18) . (a(2)) and F 2 Expressions for F 3 (a(2)) are similar to Equation (6.18) with the appropriate permutation of indices. During the limiting process, it is assumed the total number of the unit dipoles and charge carriers within the outer shell are conserved. Since generalized permittivity is directly proportional to the density of unit dipoles and unit charge carriers, an equivalent generalized surface pemfittivity e(â)V(1)(a(1)a(2)) 138 ) is related to in the limiting case as follows: = AâS(a(2)) (6.19) where S (a(2)) is the surface area of the inner ellipsoid and describes the electric nature of the surface phenomena for the homogeneous inclusion. Using Equation (6.14), it is found that for small 3(1)S((2)) g 1D(a(2)) C(a(2)) 4itC(a(2)) (620) By inserting Equations (6.12) (6.15) and (6.20) into Equation (6.7) and taldng the - limit â 0, the following approximation is obtained for the expected value of the volume integral over a homogeneous inclusion with principal semi-axis lengths (ai, a2, a3) and possessing surface phenomena described by the surface generalized pemiittivity permittivity : (E(x) - )) E(kx)dvâ) / = 4EacVinc M.j(e, a) (6.21) J pj(eb),e),Ag ,a) where Mj(c),e),,a) = - 40) + [1 A (ai, )1 S(a) V 2 - and pj(e0)cO)241)a) = ) 2 A(cci,z + (1 2 ,x 1 A(a ) ] 139 - b) 4 + 3AJ(ctl,a2) x (6.22) S(a) {F(a) + [1-Ai (ai,a )] C(a) (ai a2 a3)â} [4itC(a)]â 2 i)âą (6.23) Unlike the result obtained in Chapter Three for an inclusion without surface phenomena, it is apparent that the dielectric response of this inclusion is explicitly dependent on its dimensions. In addition, it can be seen that this approximation differs significantly from that proposed by OâKonski (1960); the effects of the surface phenomena are not expressed simply in terms of an effective inclusion permittivity. In order to test the validity of this approach, let us briefly consider two specific cases for the approximation given by Equation (6.21). First, for an ellipsoidal inclusion without surface phenomena (i.e. = 0), the result derived in Chapter Three for a homogeneous inclusion is obtained. Second, in the case of a spherical inclusion, Equation (6.21) gives J (e(x) E)EW(x)dv - () + 2r3Ea âic . (6.24) (E+2r1 2..gâ)+2e vj This result is identical to the form obtained by OâKonski (1960) and others using boundary conditions to explicitly incorporate surface phenomena. The form of this approximation can be greatly simplified by restricting the inclusion shape to a spheroid of rotation. Let a uniquely defmed by its aspect ratio, x F(a) = = al = a and b = a3; the shape of the inclusion is b/a. In this case, -[(2a2+b2)(1A(c)) 2a-2] 140 (j = 1,2), (6.25) F(a) = -[(2a2+b2)Mx) b-2] - (j = (6.26) 3), and (6.27) C(a)= (2b+?), where A((x) is the depolarization coefficient associated with the unique axis. Incorporating Equations (6.25)â(6.27) into Equation (6.21), the following estimate is obtained for the expected value of the volume integral in the case of a spheroidal inclusion: (qx) { 4fr(1) [1 -A(a)] + + - 2[1 + Ea e(o) + [1 A(c) e) + [1 A()] + + vâinc>< A(cz)] f(b,x) [1 +A(x)] e) + [1... )] (2 - - E (i)(x)clv) - + 1)(2c2 + i)â I3(b,cx)? A(x)] f(b,cL) 2A(x) 2 CL (2CL2 (6 28) + i) f3(b,a) The size dependence of the inclusion response is determined by 3(b,cC), the surface area to volume ratio for the inclusion. For oblate and prolate spheroids, 141 f is given by 13(b,c) = j (1 + x(cx i)hu12 - log [(2 + 1)1/2 + (6.29) and 13(b,a) = j-- { 1 + x(1 - 2 c,2)h/ 1 ((i sin - 2)h/2]) , (6.30) completely characterizes the effect of the respectively. Further, the quantity 13(b,cx) surface phenomena for all spheroidal particles having the same shape. Numerical Simulation To illustrate the effect of surface phenomena on the effective dielectric response of a heterogeneous media, a system containing identical spheroidal particles will be used. Utilizing Equation (6.28) to estimate the perturbation term in Equation (6.2) and employing the Lorentz approximation, the dielectric response of the resulting medium is described by 1 â { Chic X + 2E 0 4(e) - + 2[1 [1 -A(x)] )+[i +A(cL)] )+ [1 142 + A(x)] 13(b,a) ft.4a)](a2 + 1)(2cc2 + i)â f(b,c)? A(cx) + [1 A((x)] (o) + 2A(a) a 2 (2a2 - + 1)-i f(b, a) 1 , (6.31) J where cj is the volume fraction concentration of the inclusions in the system. Further, it will be assumed that the constituents and surface phenomena possess simple linear electrical rheologies of the form given in Equation (4.16) and (i) (6.30) 1(i) where and are the apparent surface permittivity and conductivity, respectively. The parameters for the system are chosen to represent a suspension of colloidal particles in a weak electrolyte (see Table 6.1). The surface phenomenon is taken to be a simple surface conduction (i.e. = 0. Farads); hence, the effects of this surface mechanism on e is determined by the quantity (b, a)j$â). Figures (6. 1)-(6.4) and Figures (6.5)-(6.8) give the computed value of i and a for this system as a function of frequency for particle shapes a = 1.000, 0.100, 0.010, and 0.001, respectively. In all cases, the fractional volume concentration of the particles was taken to be 0.0005. At this concentration level, it can be shown that identical spheroidal particles in this range of aspect ratios used can have uniformly random orientation without particle overlap. For each value of a, the response of the system is calculated for an increasing series of 13(b, a)4) values; for a constant value of 143 this represents a sequence of models where the particle size scale as expressed in terms of f(b, cx) is increasing. The presence of the surface conduction causes a significant enhancement of ic at lower frequencies for the more spherical particles (i.e. cz=1.000 and 0.100). The magnitude of the effect approaches a limiting value as f3(b, cx)i.i increases. For the more oblate particles (i.e. cx=0.010 and 0.001), the existence of the surface conduction initially causes a reduction of ic in the low frequency range that increases with increasing f3(b, ct).4â); however, this trend eventually reverses itself. In the case of cx=0.010, enhancement of K similar to that observed for the more spherical particles occurs once j3(b, is sufficiently large; conversely, enhancement of ic relative to the no surface conduction case cx)j4 used. The effect of surface conduction on a is similar in all cases; increasing f3(b, a)4 results in larger values of c is not observed for cc=0.001 for the values of f(b, over the high frequency range. Significant enhancement of a relative to the conductivity of the weak electrolyte over the lower frequencies occurs only in the case of cx=0.00 1. Conclusions A model has been presented for incorporating the effects of surface phenomena in the dielectric response of a heterogeneous medium which retains the inherent geometrical structure of this system. This is achieved by initially viewing the inclusion as a two component confocally-layered ellipsoid where the surface phenomena occur in the outer ellipsoidal shell. To obtain the approximation to a homogeneous inclusion system, the outer shell is taken to be infinitesimally thin while the total numbers of unit charge carriers and unit dipoles are conserved. The resulting estimate of the response of such an inclusion 144 shows that it is explicitly dependent on the dimensions of the inhomogeneity. This approximation is consistent with the response of an inclusion without surface phenomena and the result obtained for a spherical inclusion where surface phenomena are rigorously introduced through use of the appropriate boundary conditions. When the shape of the inclusions is restricted to spheroids of rotation, it is found that the dielectric response of the inclusion is dependent on its surface area to volume ratio. In addition, it is shown that the electric nature of the surface phenomena for the inclusions of a given shape is completely defmed by the product of this ratio and the generalized surface permittivity. Using the approximation for spheroidal inclusions, the dielectric response of a medium containing identical inclusions possessing a surface conduction has been simulated. It is found that the effect of the surface conduction changes significantly as the shape of the inclusion varies. 145 aa(in) Table 6.1: Electrolyte Inclusion 80. 5. 3 1.x10 0. Dielectric constants and conductivities used for the constituents of the model colloidal suspension. 146 80.15 80.10 cu 80.05 0. Cu 80.00 79.95 79.90 1.0 3.0 5.0 7.0 9.0 11.0 Log (Frequency (Hz)) Figure 6.1: Computed 1c as a function of frequency (plotted in log scale) for models containing identical inclusion possessing a surface conduction. Value of 3(b, c4i) used in a given model is denoted. Aspect ratio (a) of inclusion used in all models is 1.000. 147 Cu 0. Cu 79.9 1.0 3.0 5.0 7.0 9.0 11.0 Log (Frequency (Hz)) Figure 6.2: Computed ic as a function of frequency (plotted in log scale) for models containing identical inclusion possessing a surface conduction. Value of 13(b, a)I$i) used in a given model is denoted. Aspect ratio (ct) of inclusion used in all models is 0.100. 148 85 84 Cu 0. 0. Cu 83 82 81 80 79 1.0 3.0 5.0 7.0 9.0 11.0 Log (Frequency (Hz)) Figure 6.3: Computed ic as a function of frequency (plotted in log scale) for models containing identical inclusion possessing a surface conduction. Value of f3(b, a)$â) used in a given model is denoted. Aspect ratio (a) of inclusion used in all models is 0.010. 149 400 300 Cu 0. 0. Cu 200 100 0 1.0 3.0 5.0 7.0 9.0 11.0 Log (Frequency (Hz)) Figure 6.4: Computed 1c as a function of frequency (plotted in log scale) for models containing identical inclusion possessing a surface conduction. Value of f(b, a)J$1) used in a given model is denoted. Aspect ratio (cL) of inclusion used in all models is 0.001. 150 10 1 10 E 10 10 1.0 3.0 5.0 7.0 9.0 11.0 Log (Frequency (Hz)) Figure 6.5: Computed a as a function of frequency (plotted in log scale) for models containing identical inclusion possessing a surface conduction. Value of a)j4â) used in a given model is denoted. Aspect ratio (a) of inclusion 13(b, used in all models is 1.000. 151 10 I 100. 10 > / 2, 1 ioj E 0) 10 71 0.01, 0.00, 10 0.1 4 . 1.0 I âą 3.0 I 5.0 âą I 7.0 âą I âą 9.0 I 11.0 Log (Frequency (Hz)) Figure 6.6: Computed a as a function of frequency (plotted in log scale) for models containing identical inclusion possessing a surface conduction. Value of f3(b, a)J$i) used in a given model is denoted. Aspect ratio (x) of inclusion used in all models is 0.100. 152 10 I 10 10 10 1.0 3.0 5.0 7.0 9.0 11.0 Log (Frequency (Hz)) Figure 6.7: Computed a as a function of frequency (plotted in log scale) for models containing identical inclusion possessing a surface conduction. Value of 1(b, cc)4i) used in a given model is denoted. Aspect ratio (ce) of inclusion used in all models is 0.010. 153 10 I 10 10 10 1.0 3.0 5.0 7.0 9.0 11.0 Log (Frequency (Hz)) Figure 6.8: Computed a as a function of frequency (plotted in log scale) for models containing identical inclusion possessing a surface conduction. Value of a)j4â) used in a given model is denoted. Aspect ratio (x) of inclusion 1(b, used in all models is 0.001. 154 Chapter Seven Summary Microscopic phenomena can significantly affect the macroscopic physical properties of a porous rock. A mathematical model is necessary to establish and analyze the connection between a given microscopic phenomenon and the observed macroscopic behavior. Such theoretical formulations have been derived in this study for the relationship between the elastic and dielectric responses of porous rocks and certain microscopic phenomena. Contact generation has been suggested as a possible reason for the discrepancy between previous models and experimental data for the functional dependence of elastic wave velocities on confining pressure in unconsolidated materials. A model for the elastic behavior of a random sphere packing that permits the creation of new grain contacts during hydrostatic compression was derived for the first time in this thesis. Using this formulation, it was found that the closure of near-contact gaps having a relatively small mean width (i.e. on the order of 1/1000 of a sphere radius) can significantly affect the elastic behavior of the glass bead packings over the range of confining pressures employed in these experiments. Over a range of conditions that are of interest in the geophysical measurement of elastic and electromagnetic wave propagation in porous rocks, the behavior of the system is significantly affected by the geometrical configuration of the constituents on a 155 microscopic scale. Sufficient information about the microstructure of porous rocks is obtained by viewing these systems as a rock matrix into which inclusions representing the individual pore spaces are embedded. Inclusion-based formulations can be employed to determine the macroscopic physical properties of a porous rock when this description of its microstructure is adopted. Major classes of these formulations are defined by the manner in which inclusion interactions are simulated in the analysis; for each class of formulations, there is an associated microstructure that specifies the topological relationships between the various components. Pore-scale fluid distribution in partially saturated porous rocks is a microscopic phenomenon that is characterized by the geometrical configuration of the constituents within the system. Drastically different theoretical formulations were obtained for various styles of pore-scale fluid distribution while maintaining fixed constituent properties and pore structure. It was determined that the pore fluid configuration within the individual pore spaces and the pore geometries in which saturation conditions are varying are critical factors in determining the elastic and electromagnetic wave propagation in partially saturated rocks. Further, the effects of changes in saturation conditions within a particular pore geometry significantly increase as the aspect ratio of the pores decreases (i.e. as the pores become more crack-like). The mathematical expressions derived in this study were used to analyze the variation in the functional relationships observed between the level of water saturation and the physical properties of a tight gas sandstone and a gas bead packing when both of these systems were subjected to two differing saturation methods. Using simple models that consider only the basic geometrical elements that describe the pore-scale fluid 156 distribution, good agreement is obtained between the theoretical models and the experimental data. Further, it is demonstrated that a single geometrical model for the pore-scale fluid distribution can be used to simultaneously estimate the dielectric and elastic response of a partially saturated rock. However, this formulation assumes that communication between individual pores is minimal and intrapore phenomena control the physical behavior of the system. This type of formulation does not accurately predict the functional relationships observed in experimental data where these conditions are not satisfied. Surface phenomena, such as surface conduction and electrical double layers, can have a significant effect on the electrical properties of porous rocks; however, these mechanisms are explicitly excluded in the derivation of the inclusion-based formulations used to describe these systems. Further, incorporation of surface phenomena into the analysis of non-spherical inclusion shapes is very difficult. To obtain an approximation for the effects of this mechanism on the electrical response of an ellipsoidal inclusion, a model based on the limiting case of a confocally-layered inclusion was derived by assuming the conservation of unit charge carriers and unit dipoles. This approach retains the geometrical configuration inherent in the system. The behavior of the inclusion with surface phenomena is found to depend on the absolute size of the inclusion. Using this model, it is shown that the effects of the presence of a surface conduction changes significantly as the shape of the inclusion varies. The results of this study clearly demonstrate the usefulness of mathematical models in describing the effects of a given microscopic phenomenon on the macroscopic physical 157 properties of a porous rock. It is significant to note that such an analysis does not always require a detailed description of the porous rock system; good agreement between the observed behavior and a theoretical formulation can be obtained by utilizing a relatively simple model. However, much basic work still remains to be done in this field. Not only do models for other microscopic phenomena and macroscopic physical properties need to be developed, the use of these formulations to rigorously analyze experimental measurements in the context of inverse problem theory must be further explored. Therefore, I foresee much significant work occurring in this field over both the near-term and distant future. 158 References Avellaneda, M., 1987, âIterated Homogenization, Differential Effective Medium Theory and Applications,â Commun. Pure. Appl. Math., v. 40, pp. 527-554. Batchelor, G. K., and OâBrien, R. W., 1977, âThermal or Electrical Conduction through a Granular Material,â Proc. R. Soc. Lond., v. A355, pp. 313-333. Bernal, J. D., and Mason, J., 1960, âCo-ordination of Randomly Packed Spheres,â Nature, v. 188, pp. 910-911. Biot, M. A., 1956, âTheory of Propagation of Elastic Waves in a Fluid-Saturated Porous Solid. I. Low Frequency Range,â 3. Acoust. Soc. Am., v. 28, pp. 168-178. BĂ¶ttcher, C. 3. F., 1942, âZur Theorie der Inneren Elektrischen FeldstĂ€rke,â Physica, v. 9, pp. 937-944. BĂ¶ttcher, C. J. F., 1945, âThe Dielectric Constant of Crystalline Powders,â Rec. Tray. Chim. Pays-Bas, V. 64, pp. 47-51. Bourbie, T., and Zinszner, B., 1984, âSaturation Methods and Attenuation versus Saturation Relationships in Fontainebleau Sandstone,â Presented 54 th Ann. Intemat. Mtg., Soc. Explor. Geophys., Atlanta, Expanded Abstract, pp. 344-347. Brace, W. F., Silver, E., Hadley, K., and Goetze, C., 1972, âCracks and Pores: A Closer Look,â Science, v. 178, pp. 162-169. Brandt, H., 1955, âA Study of the Speed of Sound in a Porous Granular Media,â 3. Appi. Mech., v. 22, pp. 479-486. Bruggeman, D. A. G., 1935, âBerechnung Verschiedener Physikalischer Konstanten von Heterogenen Substanzen, I. DielekirizitĂ€tskonstanten und LeitfĂ€higkeiten der MischkĂ¶rper aus Isotropen Substanzen,â Annalen der Physik 5, vol. 24, pp. 636-664. Burger, H. C., 1919, âDas Leitvermogen VerdĂŒnnter Mischkristallfreier Legierungen,â Physik. Zeitschr., v. 20, pp. 73-75. Burington, R. S., 1965, Handbook of Mathematical Tables and Formulas, McGraw-Hill, 159 New York. Cheng, C. H., 1978, âSeismic Velocities in Porous Rocks: Direct and Inverse Problems,â Ph.D. dissertation, Massachusetts Institute of Technology, Cambridge, Mass. Cheng, C. H., and ToksĂ¶z, N. M., 1979, âInversion of Seismic Velocities for the Pore Aspect Ratio Spacetrum of a Rock,â 3. Geophys. Res., v. 84, pp. 7533-7543. Clark, V. A., Tittman, B. R., and Spencer, T. W., 1980, âEffect of Volatiles on Attenuation (Q) and Velocity in Sedimentary Rocks,â J. Geophys. Res., v. 85, pp. 5 1905198. Clausius, R. 3. E., 1879, Die Mechanische WĂ€rmetheorie, v. 2, Die Mechanische Behandlung der ElectricitĂ€t, Vieweg, Braunschweig. Cohen, M. H., 1981, âScale Invariance of the Low-Frequency Electrical Properties of Inhomogeneous Materials,â Geophysics, v. 46, pp. 1057-1059. De La Rue, R. E., and Tobias, C. W., 1959, âOn the Conductivity of Dispersions,â 3. Electrochem. Soc., v. 106, pp. 827-833. de Loor, G. P., 1953, âMethod of Obtaining Information on the Internal Dielectric Constant of Mixtures,â Appl. Sci. Res. B., v. 3, pp. 479-482. de Loor, G. P., 1964, âDielectric Properties of Heterogeneous Mixtures with a Polar Constituent,â Appl. Sci. Res. B., v. 11, pp. 310-320. de Loor, G. P., 1968, âDielectric Properties of Heterogeneous Mixtures Containing Water,â J. Microwave Power, v. 3, pp. 67-73. Deresiewicz, H., 1958, âStress-Strain Relations for a Simple Model of a Granular Medium,â 3. Appl. Mech., v. 25, pp. 402-405. Digby, P. 3., 1981, âThe Effective Elastic Moduli of Porous Granular Rocks,â 3. Appl. Mech., v. 48, pp. 803-808. Domenico, S .N., 1976, âEffect of Brine-Gas Mixture on Velocity in an Unconsolidated Sand Reservoirs,â Geophysics, v. 41, pp. 882-894. 160 Domenico, S. N., 1977, âElastic Properties of Unconsolidated Porous Sand Reservoirs,â Geophysics, v. 42, pp. 1339-1368. Duffy, 3., and Mindlin, R. D., 1957, âStress-Strain Relations of a Granular Medium,â 3. Appi. Mech., v. 24, pp. 585-593. Dukhin, S. 5., 1971, âDielectric Properties of Disperse Systems,â in Surface and Colloid Science, v. 3, ed. E. Matijevic, Wiley, New York, pp. 83-165. Dukhin, S. S., and Shilov, V. N., 1974, Dielectric Phenomena and the Double Layer in Disperse Systems and Polyelectrolytes, Wiley, New York. Fixman, M., 1980, âCharge Macromolecules in External Fields. I. The Sphere,â 3. Chem. Phys., v. 72, pp. 5177-5186. Fixman, M., 1983, âThin Double Layer Approximation for Electrophoresis and Dielectric Response,â J. Chem. Phys., v. 78, pp. 1483-1491. Foster, A. G., 1932, âThe Sorption of Condensible Vapors by Porous Solids. Part I. The Applicability of the Capillary Theory,â Trans. Faraday Soc., v. 28, pp. 645-657. Fricke, H., 1924, âA Mathematical Treatment of the Electric Conductivity and Capacity of Disperse Systems: I. The Electric Conductivity of a Suspension of Homogeneous Spheroids,â Phys. Rev., v. 24, pp. 575-587. Fricke, H., 1953, âThe Maxwell-Wagner Dispersion in a Suspension of Ellipsoids,â J. Phys. Chem., v. 57, pp. 934-937. Fuller, B. D., and Ward, S. H., 1970, âLinear System Description of the Electrical Parameters of a Rock,â I.E.E.E. Trans. Geosci. Elect., v. GE-8, pp. 7-18. Gassmann, F., 1951a, âElastic Waves through a Packing of Spheres,â Geophysics, v. 16, pp. 673-685. Gassmann, F., 195 ib, âUeber die ElastizitĂ€t PorĂ¶ser Medien,â Vierteljahrsschr. Naturforsch. Ges. Zurich, v. 96, pp. 1-2 1. Gradshteyn, I. S., and Ryzhik, I. M., 1965, Table ofIntegrals, Series, and Products, Academic Press, New York. 161 Gregory, A. R., 1976, âFluid Saturation Effects on Dynamic Elastic Properties of Sedimentary Rocks,â Geophysics, v. 41, pp. 895-921. Haines, W. B., 1930, âStudies in the Physical Properties of Soils. V. The Hysteresis Effect in Capillary Properties and the Modes of Moisture Distribution Associated Therewith,â J. Agric. Sci., v. 20, pp. 97-116. Knight, R., and Endres, A., 1990, âA New Concept in Modeling the Dielectric Response of Sandstones: Defining a Wetted Rock And Bulk Water System,â Geophysics, v. 55, pp. 586-594. Knight, R. and Nolen-Hoeksema, R., 1990, âA Laboratory Study of the Dependence of Elastic Wave Velocities on Pore Scale Fluid Distribution,â Geophy. Res. Lett., v. 17, pp. 1529-1532. Knight, R. and Nur, A., 1987a, âGeometrical Effects in the Dielectric Response of Partially Saturated Sandstones,â Log Anal., v. 28, pp. 513-519. Knight, R. and Nur, A., 1987b, âThe Dielectric Constant of Sandstones, 60 kHz to 4 MHz,â Geophysics, v. 52, pp. 644-654. Kuster, G. T. and ToksĂ¶z, M. N., 1974, âVelocity and Attenuation of Seismic Waves in Two-Phase Media: Part I. Theoretical Formulations,â Geophysics, v. 39, pp. 587-606. Landau, L. D. and Lifshitz, E. M., 1960, Electrodynamics of Continuous Media, Pergamon, Oxford. Landauer, R., 1952, âThe Electric Resistivity of Binary Metallic Mixtures,â J. Appl. Phys., v.23, pp.779-784. Landauer, R., 1978, âElectric Conductivity in Inhomogeneous Media,â in Electrical Transport and Optical Properties in Inhomogeneous Media, eds. J. C. Garland and D. B. Tunner, Am. Inst. Phys., New York, pp. 2-43. Longeron, D. F., Arguad, M. J., and Ferand, 3. P., 1989, âEffect of Overburden Pressure, Nature, and Microscopic Distribution of the Fluids on Electrical Properties of Rock Samples,â SPE Form. Eval., v. 4, pp. 194-202. Lorentz, H. A., 1880, âUeber die Beziehung Zwischen der Fortpflanzungsgeschwindigkeit 162 des Lichtes der Korperdichte,â Annalen der Physik, v. 9, pp. 641-665. Lorenz, L. V., 1880, âUeber die Refractionsconstante,â Annalen der Physik, v. 11, pp. 70-103. Meredith, R. E., and Tobias, C. W., 1961, âConductivities in Emulsions,â J. Electrochem. Soc., v. 108, pp. 286-290. Meredith, R. E., and Tobias, C. W., 1962, âConduction in Heterogeneous Systems,â in Advances in Electrochemistry and Electrochemcical Engineering, v.2, ed. C. W. Tobias, Interscience, New York, pp. 15-47. Miles, J. B., and Robertson, H. P., 1932, âThe Dielectric Behavior of Colloidal Particles with an Electrical Double Layer,â Phys Rev., v. 40, pp. 583-591. Milton, G. W., 1984, âCorrelation of the Electromagnetic and Elastic Properties of Compsites and Microgeometries Corresponding with Effective Medium Approximations,â in Physics and Chemistry of Porous Media, eds. D. L. Johnson and P. N. Sen, Am. Inst. Phys., New York, pp. 66-77. Milton, G. W., 1985, âThe Coherent Potential Approximation is a Realizable Effective Medium Scheme,â Commun. Math. Phys., v. 99, pp. 463-500. Mossotti, 0. F., 1850, âDiscussione Analitica Sullâinfluenza che Lâazione di un Mezzo Dieletthco ha Sulla Distribuzione DellâelettricitĂ alla Superfizie di piĂŒ Corpi Electtrici Disseminati in Esso,â Mem. di Matem. di Fisica Soc. hal. Sci. Modena 2, v. 24, pp. 4974. Murphy, W. F., 1982a, âEffects of Microstructure and Pore Fluids on Granular Sedimentary Materials,â Ph.D. dissertation, Stanford University, Palo Alto, Calif. Murphy, W. F., 1982b, âEffects of Partial Water Saturation on Attenuation in Massilon Sandstone and Vycor Porous Glass,â J. Acoust Soc Am., v. 71, pp. 1458-1468. Murphy, W. F., 1984, âAcoustic Measures of Partial Gas Saturation in Tight Gas Sandstones,â J. Geophys. Res., v. 89, pp. 11549-11559. Murphy, W. F., 1985, âSonic and Ultrasonic Velocities: Theory versus Experiment,â Geophys. Res. Lett., v. 12, pp. 85-88. 163 Norris, A. N., Callegari, A. J., and Sheng, P., 1985a, âA Generalized Differential Effective Medium Theory,â J. Mech. Phys. Solids., v. 33, pp. 525-543. Norris, A. N., Sheng, P., and Callegari, A. J., 1985b, âEffective-Medium Theories for Two-Phase Dielectric Media,â J. Appi. Phys., v. 57, pp. 1990-1996. Nur, A., and Simmons, G., 1969, âThe Effects of Saturation on Velocity in Low Porosity Rocks,â Earth Plan. Sci. Lett., v. 7, pp. 183-193. OâKonski, C. T., 1960, âElectric Properties of Macromolecules. V. Theory of Ionic Polarization in Polyelectrolytes,â 3. Phys. Chem., v. 64, pp. 605-619. Onsager, L., 1936, âElectric Moments of Molecules in Liquids,â 3. Am. Chem. Soc., v. 58, pp. 1486-1493. Papoulis, A., 1965, Probability, Randt,m Variables, and Stochastic Processes, McGrawHill, New York. Pearce, C. A. R., 1955, âThe Electhcal Conductivity and Permittivity of Mixtures, with Special Reference to Emulsions of Water in Fuel Oil,â Brit. I. Apple. Phys., v. 6, pp. 113-120. Polder, D., and Van Santen, 3. H., 1946, âThe Effective Permeability of Mixtures and Solids,â Physica, v. 12, pp. 257-271. Poley, 3. Ph., Nootebom, J. J., and de Waal, P. 3., 1978, âUse of V.H.F. Dielectric Measurements for Borehole Formation Analysis,â Log Analyst, v. 20, n. 3, pp. 8-30. Schurr, J. M., 1964, âOn the Theory of the Dielectric Dispersion of Colloidal Particales in Electrolyte Solution,â J. Phys. Chem., v. 68, pp. 2407-2413. Schwarz, G., 1962, âA Theory of the Low-Frequency Dielectric Dispersion of Colloidal Particles in Electrolyte Solution,â 3. Phys. Chem., v.66, pp. 2636-2642. Sihvola, A. H., and Kong, J. A., 1988, âEffective Permittivity of Dielectric Mixtures,â IEEE Trans. Geosci. Remote Sen. v. 26, pp. 420-429; Erratum, 1989, v. 27, pp. 101102. , 164 Sprunt, E. S., and Brace, W. F., 1974, âDirect Observations of Microcavities in Crystalline Rocks,â Tnt. J. Rock Mech. Mi Sci. and Geomech., Abst., v. 11, pp. 139150. Stepin, L. D., 1965, âDielectric Permeability of a Medium with Non-Uniform Ellipsoidal Inclusions,â Soy. Phys.-Tech. Phys., v. 10, pp. 768-772. Stratton, J. A., 1941, Electromagnetic Theory, McGraw-Hill, New York. Stroud, D., 1975, âGeneralized Effective-Medium Approach to the Conductivity of an Inhomogeneous Material,â Phys. Rev. B, v. 12, pp. 3368-3373. Swanson, B. F., 1979, âVisualizing Pores and Non-wetting Phase in Porous Rock,â J. Pet. Tech., v. 31, pp. 10-18. Taylor, L. S., 1965, âDiectric Properties of Mixtures,â I.E.E.E. Trans. Anten. Propag., v. AP-13, pp. 943-947. Timur, A., Hempkins, W. B., and Weinbrandt, R. M., 1971, âScanning Electron Microscopic Study of Pore Systems in Rocks,â J. Geophys. Res., v. 76, pp. 4932-4948. ToksĂ¶z, M. N., Cheng, C. H., and Timur, A., 1976, âVelocities of Seismic Waves in Porous Rocks,â Geophysics, v. 41, pp. 621-645. van Beek, L. K. H., 1967, âDielectric Behaviour of Heterogeneous Systems,â in Progress in Dielectrics, v. 7, ed. J.B.Birks, Heywood, London, pp. 69-114. Wagner, K. W., 1914, âErkiarung der Dielektrischen Nachwirkungen auf Grund Maxwellscher Vorstellungen,â Arch. Elektrotech., v. 2, pp. 37 1-387. Walls, J. D., 1982, âTight Gas Sands-Permeability, Pore Structure, and Clay,â J. Pet. Tech., v. 34, pp. 2708-2714. Walton, K., 1975, âThe Effective Elastic Moduli of Model Sediments,â Geophys. J. R. Astro. Soc., v. 43, pp. 293-306. Walton, K., 1978, âThe Oblique Compression of Two Elastic Spheres,â J. Mech. Phys. Solids, v. 26, pp. 139-150. Walton, K., 1987, âThe Effective Elastic Moduli of a Random Packing of Spheres,â 3. 165 Mech. Phys. Solids, v. 35, pp. 2 13-226. White, J. E., 1965, Seismic Waves: Radiation, Transmission and Attenuation, McGrawHill, New York. Wyllie, M. R. J., Gregory, A. R., and Gardner, L. W., 1956, âElastic Wave Velocities in Heterogeneous and Porous Media,â Geophysics, v. 21, pp. 4 1-70. Yonezawa, F., and Cohen, M. H., 1983, âGranular Effective Medium Approximation,â J. Appl. Phys., v. 54, pp. 2895-2899. 166 Appendix One Constitutive Relationships and Generalized Permittivity The electromagnetic behavior of a homogeneous, isotropic medium is generally expressed in terms of its dielectric permittivity magnetic permeability (pP) (E), electrical conductivity (a), and which describe the capacitive, conductive, and inductive properties, respectively, of the substance. For materials exhibiting a linear electromagnetic behavior, these quantities are defined by the following relationships for the time harmonic fields in complex notation: D =cE, (Au) JfGE, (A12) BtH, (AL3) and where B , D , E, H, andJ are the magnetic induction, electric displacement, electric field strength, magnetic field strength, and electric free current density, respectively. The quantities c, a , and t are the Fourier transforms of the impulse response for each of the 167 respective properties. In general, e, a, and jt are complex functions of the form = eâ i eâ. (A 1.4) iaâc (AL5) .tâ. (A L6) - 0=0â- and = iiâ - where the superscriptsâ andâ denote the real and imaginary parts. Assuming causality, the real and imaginary parts are related through the Hilbert transform. The electrical properties are generally measured by observing the electric current density which results from an applied electric field. Assuming that the effects of magnetic induction may be neglected, this observed current density is the total current density Jt where J Jd+Jf (At7) and JdiD 168 (A18) is the displacement current density associated with the capacitive behavior of the medium (Fuller and Ward, 1970). Therefore, it is useful to revise the constitutive relationship given by Equation (A 1.2) to the form J= aRE. (A1.9) where (A1.lo) is the generalized electrical conductivity. Conversely, the generalized dielectric permittivity 8 g can be defined in term of the generalized electric displacement D as follows: DEgE (A1.1l) where D =D + Jf 1 (io) and 169 (Au2) Eg= O. 1 EioY (AL13) A complete description of the electrical properties of a linear material can be accomplished by employing either generalized constitutive relationship given by Equations (A 1.9) and (A 1.11). The real part of e and (Yg are referred to as the apparent dielectric permittivity Ea and the apparent electrical conductivity Ga, respectively. It is these quantities that are generally reported as the experimentally measured dielectric permittivity and electhcal conductivity. Further, coupling of the capacitive and conductive properties occurs at the discrete interfaces within a heterogeneous medium. Hence, generalized constitutive relationship provide an effective means for incorporating both types of phenomena into the inclusion-based formulation for the electrical properties of porous rocks. 170 Appendix Two Evaluation of the Depolarization Coefficient Associated with the Unique Principal Axis of the Inner Spheroid within a Confocally-Layered Inclusion Evaluation of the expression describing the dielectric response of a confoc ally-layered inclusion requires the determination of B ((i), x), the depolarization coefficient associated with the unique axis of the inner spheroid. For a non-degenerate spheroidal pore shape (i.e. 0< a<1), this is accomplished by determining the aspect ratio of the inner spheroid (f). Let (ai,a ) and 2 2 ,b be the major and minor semi-axis length pair for the pore 1 (b ) space and inner spheroid, respectively; then, 1 a (A2.1) and (A2.2) Since the two spheroids are confocal, then the semi-axis lengths are related by -(bi) 2 (ai) = 2 -(b (a . ) 171 (A2.3) Using Equations (A2. 1) and (A2.2), Equation (A2.3) can be rewritten as i2 2 1-a J 1 b (A2.4) The fractional volume of the non-wetting pore fluid phase within the inclusion is given by 2 (b b 2 ) 1 (A2.5) 2 (ai) a 2 Again using Equations (A2.l) and (A2.2), Equation (A2.5) can be expressed as p2/3 = 213 (ai)2 (ccs(1)) (A2.6) Combining Equations (A2.4) and (A2.6), one obtains (1))m (a 5 where y= f (A2.7) 13 Equation (A2.7) is the normal form of a cubic system which is discussed in Burington (1965). For a non-degenerate oblate spheroidal pore, this equation possesses only one real-valued root that is given by = 2 1 3 [2()1â2cot(2Ă)] , â 172 (A2.8)) where = tan1[(tan (p)113], (A2.9) (A2. 10) = and = 2 1-c (A2.11) (c s(1))213 Once j3 is determined, B(s(â),ĂĄ)is calculated using Equations (4.1) and (4.3). In the degenerate cases of the sphere and the flat disc, respectively. 173 B(s(1),X) is taken to be and 1,
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A theoretical treatment of microscopic phenomena in...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A theoretical treatment of microscopic phenomena in porous rocks Endres, Anthony Lee 1992
pdf
Page Metadata
Item Metadata
Title | A theoretical treatment of microscopic phenomena in porous rocks |
Creator |
Endres, Anthony Lee |
Date Issued | 1992 |
Description | This thesis makes five significant contributions to the theoretical analysis of the connection between microscopic phenomena and the macroscopic physical properties of porous rocks. A random sphere packing model permitting contact generation during hydrostatic compression is derived. It is demonstrated that the closure of near-contact gaps with an extremely small mean width significantly alters the elastic properties of granular media from those predicted by previous models in the pressure range used in laboratory measurements. Generalized forms of the inclusion-based formulations are obtained; major classes of these formulations are defined according to the manner in which interactions between inclusions in the heterogeneous system are simulated. Each class possesses an associated microstructure that determines the topological relationship between the various components. Inclusion-based formulations are obtained that describe the effects of pore-scale fluid distribution on the dielectric and elastic responses of a partially saturated rock. It is found that the pore fluid configuration within the individual pores and the pore geometries in which saturation conditions are varying are critical factors in determining the dielectric and elastic properties of partially saturated rocks. In addition, it is observed that the effect of saturation condition variations in a particular pore geometry increases as the pore shape becomes more crack-like. The theoretical formulations describing the effects of pore-scale fluid distribution are used to analyze experimental data for a partially saturated tight gas sandstone and glass bead packing. Simple models that incorporate only the basic geometrical elements of the resulting pore-scale fluid distributions accurately predict the experimental data. It is also found that the same geometrical model can be used simultaneously to estimate the dielectric and elastic responses of a partially saturated porous medium. The effect of surface phenomena (e.g. electrical double layers and surface conduction) at the solid-fluid interface are incorporated into inclusion-based formulations of the dielectric response by employing the limiting case of a confocally-layered ellipsoid. It is found that the effect of surface phenomena varies as a function of both inclusion size and shape. |
Extent | 2247220 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2008-12-23 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0052955 |
URI | http://hdl.handle.net/2429/3278 |
Degree |
Doctor of Philosophy - PhD |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1992-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_1992_spring_endres_anthony_lee.pdf [ 2.14MB ]
- Metadata
- JSON: 831-1.0052955.json
- JSON-LD: 831-1.0052955-ld.json
- RDF/XML (Pretty): 831-1.0052955-rdf.xml
- RDF/JSON: 831-1.0052955-rdf.json
- Turtle: 831-1.0052955-turtle.txt
- N-Triples: 831-1.0052955-rdf-ntriples.txt
- Original Record: 831-1.0052955-source.json
- Full Text
- 831-1.0052955-fulltext.txt
- Citation
- 831-1.0052955.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0052955/manifest