Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

The characterization of porous rocks with nuclear magnetic resonance and the confocal laser scanning… Chapman, Alice Elizabeth 1992

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_1992_spring_chapman_alice.pdf [ 6.11MB ]
Metadata
JSON: 831-1.0052496.json
JSON-LD: 831-1.0052496-ld.json
RDF/XML (Pretty): 831-1.0052496-rdf.xml
RDF/JSON: 831-1.0052496-rdf.json
Turtle: 831-1.0052496-turtle.txt
N-Triples: 831-1.0052496-rdf-ntriples.txt
Original Record: 831-1.0052496-source.json
Full Text
831-1.0052496-fulltext.txt
Citation
831-1.0052496.ris

Full Text

THE CHARACTERIZATION OF POROUS ROCKS WITH NUCLEARMAGNETIC RESONANCE AND THE CONFOCAL LASER SCANNINGMICROSCOPEbyALICE ELIZABETH CHAPMANB.Sc., The University of British Columbia, 1989A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinTHE FACULTY OF GRADUATE STUDIES(Department of Geological Sciences)We accept this thesis as conformingto the required standardTHE UNIVERSITY OF BRITISH COLUMBIASeptember, 1992© Alice Elizabeth Chapman, 1991In presenting this thesis in partial fulfilment of the requirements for an advanceddegree at the University of British Columbia, I agree that the Library shall make itfreely available for reference and study. I further agree that permission for extensivecopying of this thesis for scholarly purposes may be granted by the head of mydepartment or by his or her representatives. It is understood that copying orpublication of this thesis for financial gain shall not be allowed without my writtenpermission.Department of ;77-L() (-7/C JfrLThe University of British ColumbiaVancouver, CanadaDate Cc /3,, /9gQDE-6 (2/88)ABSTRACTThe characterization of porous media with the confocal laser scanning microscope (CLSM)and nuclear magnetic resonance (NMR) was investigated.Pore casts (epoxy replicas) of one sandstone, Berea 100, and two carbonates, S55 andS77, were imaged with the CLSM. Isolated two-dimensional images and serial sections, a seriesof two-dimensional images separated by a specific vertical distance, were collected. The isolatedimages were analysed and mean pore radii of 25.3, 17.5, and 90.2p.m were determined for Berea100, S55 and S77 respectively.The serial sections were reconstructed to form three-dimensional images of the porenetworks. However, the pores in the samples were large compared to the depth of the pore casts,resulting in images with insufficient depth to allow three-dimensional pore size analysis.Spin-lattice (T1)measurements of twenty-two saturated carbonates (samples S31 to S49,S54, S55, S77) and six partially saturated carbonates were analysed to determine pore sizes andfluid distributions respectively. Pore sizes were calculated using the “two-fractions in fastexchange” model and the slow diffusion model of Brownstein and Tarr (1979). Samples S31 toS49, with pore radii of less than 41.Lm, were within the fast diffusion regime and thus the “two-fractions in fast exchange model” produced reliable pore sizes. Samples S54, S55 and S77, withmean pore radii of 5.1, 19.3, and 85.3.tm, were not within the fast diffusion regime and so theaccurate analysis of data for these samples required the slow diffusion model.The spin-lattice time constants calculated for the partially saturated samples were assumedto be proportional to the sizes of the saturated pores. Thus it was concluded that the larger poreswithin the samples were drained preferentially; however, large pores which were accessiblethrough smaller pore throats remained saturated to low saturations.Further investigation into these methods could focus on the production of pore casts toenable reconstruction of an image with sufficient depth for three-dimensional analysis or on thedevelopment of a pore model which more closely approximate the shape of the pores.IITABLE OF CONTENTSABSTRACT .iiTABLE OF CONTENTS mLIST OF TABLES viLIST OF FIGURES viiACKNOWLEDGEMENTS xiiCHAPTER ONE Introduction 1CHAPTER TWO Characterization of Porous Rocks with theConfocal Laser Scanning Microscope 4INTRODUCrION 4Characterizing Porous Media 4Mercury Injection 5Gas AdsorptionfDesorption 5Electrical and Optical Microscopy 6THE CONFOCAL LASER SCANNING MICROSCOPE (CLSM) 9EXPERIMENTAL PROCEDURES 11Sample Preparation 11Materials 12Pore Cast Production 14Impregnation Techniques 15Wet Impregnation 15Vacuum ImpregnationImmersion Technique 16Vacuum ImpregnationLowering Technique 16Etching 17Double Pore Casts 18DATA COLLECfION AND ANALYSIS 18RESULTS 24DISCUSSION AND CONCLUSIONS 48I I ICHAPTER THREE Characterizing Porous Rocks withNuclear Magnetic Resonance 54INTRODUCTION 54HISTORY OF NMR IN PETROLEUM GEOLOGY 55Nuclear Magnetic Logging (NML) 55Relationships Between NMR and Physical Properties 56Pore Sizes 56Fluid Distributions 57THEORY 57Longitudinal Relaxation Time (Ti) 60Transverse Relaxation Time (T2) 61The NMR Response of Fluids Confmed in Rocks 62Slow Diffusion 63Fast Diffusion 65EXPERIMENTAL PROCEDURES 67Sample Preparation 67Sample Characterization 68Porosity and Permeability 68Mineralogy 69Samples S31 to S49 72Samples S54,S55, andS77 72DATA COLLECHON 76Sample Preparation 76Fully Saturated Samples 76Partially Saturated Samples 77NMR Measurements 78Spin-Lattice Decay Curves 78Spin-Spin Decay Curves 79DATA ANALYSIS 79Data Inversion 80Spin-Lattice Relaxation Time (Ti) Distributions 80Pore Sizes - Slow Diffusion Regime 82Porosity Determination 83Pore Sizes - Fast Diffusion Regime 84RESULTS 84Fully Saturated Samples 84Ti Distributions 84Porosity 85Pore Sizes 90Fast Diffusion 90Slow Diffusion 92ivPartially Saturated Samples .102Fluid Distributions 102Fluid Saturations 102DISCUSSION AND CONCLUSIONS 108CHAPTER FOUR Summary 113INTRODUCTION 113CONFOCAL MICROSCOPY 113NMR 116Determination of Pore Sizes 116Fluid Distributions 119CONCLUSIONS 120REFERENCES 121APPENDIX 1 BioRad MRC-500 Image File Format 127APPENDIX 2 Image Analysis Programs 129APPENDIX 3 Berea 100 Confocal Pore Size Distributions 165APPENDIX 4 S55 Confocal Pore Size Distributions 174APPENDIX 5 S77 Confocal Pore Size Distributions 185APPENDIX 6 Carbonate Samples NMR Ti and Pore Size Distributions 193APPENDIX 7 NMR Ti distributions and M(0) vs. Sw graphs forPartially Saturated carbonate samples 238VLIST OF TABLES2.1 Spurr Low-Viscosity Embedding Media 132.2 Berea 100 Mean Pore Sizes 472.3 Carbonate (S55) Mean Pore Sizes 472.4 Carbonate (S77) Mean Pore Sizes 482.5 Mean Pore Sizes 503.1 Porosity and Permeability 703.2 Stains for Differentiating Between Calcite and Dolomite 71viLIST OF FIGURES2.1.20Pixel intensities in rows 40 and 250 of image ABR1. The intensity value 50 isindicated as possible separation value between pore pixels and solid pixels2.2a 21Binary versions of Berea 100 image ABR1 with pore/solid separation value of 40.Black = solid, and white = pore space.2.2b 22Binary versions of Berea 100 image ABR1 with pore/solid separation value of 50.Black = solid, and white pore space.2.2c 23Binary versions of Berea 100 image ABR1 with pore/solid separation value of 60.Black = solid, and white = pore space.2.3 25Original version of a confocal laser scanning microscope image of Berea 100. Theimage is called ABR1 and the pore space corresponds to the brightest areas of theimage.2.4a 26CLSM image of vugs in the carbonate sample, S55. The image is called F78 andthe brightest areas of the image represent porosity.2.4b 27CLSM image of vugs in the carbonate sample, S55. The image is called T78 andthe brightest areas of the image represent porosity.2.5 29Original version of a confocal laser scanning microscope image of carbonate sampleS77. The image is called Fl and the pore space corresponds to the brightest areas ofthe image.2.6a 31Binary version of carbonate S55 CLSM image F78 (Figure 2.4a) produced using asolid/pore separation value of 50. White = porosity, and Black = solid.2.6b 32Binary version of carbonate S55 CLSM image F78 (Figure 2.4a) produced using asolid/pore separation value of 100. White = porosity, and Black = solid.2.7 33Binary version of carbonate S77 CLSM image Fl (Figure 2.5) produced using asolid/pore separation value of 15. White = porosity, and Black = solid.vii2.8.35Binary versions of Berea 100 image ABR1 with pore throats separated from thepores and the centers from which the sizes of the pores and pore throats will bemeasured indicated by a contrasting colour within the void space. Black = solid,the White = porosity, and the dotted areas = pore throats.2.9a 37Examples of the Berea 100 pore size distributions. The upper plot represents thepore size distribution of image ABR 1, the lower of image HBR 1.2.9b 38Berea 100 cumulative pore size distributions. The upper plot represents thesummation of the distributions of 5 images, the lower of 16 images.2.lOa 39Examples of the carbonate (S55) pore size distributions. The upper plot representsthe pore size distribution of image F78, the lower of image N78.2.lOb 40Carbonate (S55) cumulative pore size distributions. The upper plot represents thesummation of the distributions of 5 images, the lower of 20 images.2.lla 42Examples of the carbonate (S77) pore size distributions. The upper plot representsthe pore size distribution of image Gi, the lower of image Ji.2.llb 43Carbonate (S77) cumulative pore size distributions. The upper plot represents thesummation of the distributions of 5 images, the lower of 13 images.2.12 44Porosity (the ratio of pore pixels to the total number of pixels) of the Berea 100images.2.13 45Porosity (the ratio of pore pixels to the total number of pixels) of the carbonate S55images. The upper distribution represents the binary images produced using a lowintensity separation value, and the lower distribution represents those imagesproduced with a high intensity separation value.2.14 46Porosity (the ratio of pore pixels to the total number of pixels) of the carbonate S77images.2.15 51Berea 300 capillary pressure pore size distributions. (from Crocker et al., 1983viii2.16.53Consecutive optical sections of Berea 100. The spacing was 9.6 microns betweenthe sections. From right to left, the images are the second, third and final (fourth)images collected. ABR1 was the first of this set (Figure 2.4). The bright areas inthe images are the pores, and the darker areas represent the solid.3.1 59Nuclear Magnetic Resonance Reference Frame3.2 73Partial dolomitization in Elkton Member samples, S31 to S49. Dolomite is blue andusually rhombic, calcite is red or pink, and extremely small crystals consisting ofcarbonate and silicate minerals appear as dark areas.3.3 74Examples of crinoid stems and ossicles in Elkton Member samples, S31 to S49.Dolomite is blue and usually rhombic, calcite is red or pink, and extremely smallcrystals consisting of carbonate and silicate minerals appear as dark areas. Lowerphoto contains a cross section of a crinoid stem3.4 75Examples of bryozoan and brachiopods which occur in Elkton Member samples,S31 to S49. Dolomite is blue and usually rhombic, calcite is red or pink, andextremely small crystals consisting of carbonate and silicate minerals appear as darkareas.3.5 86Sample S39 Ti distributions.(a) discrete and (b) continuous.3.6 87Sample S31 Ti spectra.(a) discrete and (b) continuous.3.7 88Sample S77 Ti spectra. (a) discrete and (b) continuous, interpolate a coefficient of2.25x10-9 m2/s at 40°C, the temperature of the sample chamber in the NMRspectrometer used to collect the data for this study.3.8 89Correlation of the integral over the the discrete and continuous Ti distributions(M(0)) with measured porosity.3.9 91100,000 ppm NaCl brine relaxation time.3.10 93(a) Sample S39 Ti and fast diffusion pore size distributions. (b) S3 i fast diffusionpore size distribution.3.11 94(a) Sample S77 Ti and fast diffusion pore size distributions. (b) S55 fast diffusionpore size distribution.ix3.12.95Slow diffusion method results of inverting theoretical data produced using atheoretical pore size spectrum similar to that of S77. Variation in the pore sizes is aresult of using two time ranges to calculate the theoretical decay curve. The uppergraph was produced with data over the same time range as the data for S77 wascollected, the lower with an extended range.3.13 96Slow diffusion method results of inverting theoretical data based on a theoreticalpore size spectrum similar to that of S77. The pore sizes resulting when 1 and 5%noise are added to the theoretical data are shown in the upper and lower graphsrespectively.3.14 99Slow diffusion pore size distributions. (a) S43 (b) S313.15 100Slow diffusion pore size distributions. (a) S55 (b) S773.16 101Examples of the variation in the pore sizes determined using the slow diffusionmethod when the value of the bulk fluid relaxation time, Tib, is varied.3.17 103Examples of the variation in the pore sizes determined using the slow diffusionmethod when the value of the self diffusion coefficent of the bulk fluid, D, isvaried.3.18 104Examples of the variation in the pore sizes determined using the slow diffusionmethod when the value of the surface sink parameter, p, is varied.3.19 105Example of the decrease in amplitude and relaxation time with saturation forpartially saturated samples.3.20 106Example of the increase of the amplitude of higher relaxation times at lowsaturations.3.21 107Examples of the correlation between fluid saturation and the total integral over theTi distribution (M(0)).3.22 iiiComparisons of the pore size spectra calculated from the NMR data with the fastand slow diffusion methods with the spectra determined using the confocal laserscanning microscope.x4.1.114Berea 100 pore size distributions. The upper distribution is a capillary pressure poresize distribution (after Crocker et aL, 1983) and the lower is the distributiondetermined from 16 confocal laser scanning microscope images.4.2 115Carbonate (S55 and S77) cumulative pore size distributions. The upper plotrepresents the summation of the distributions of 20 S55 images, the lower of 13S77 images.4.3 118Ti and fast diffusion pore size disthbutions. (a) Sample S39 (b) Sample S31xiACKNOWLEDGEMENTSI would like to express my appreciation of the enthusiasm and support provided by mysupervisor, Dr. Rosemary Knight. I would also like to thank Michael Knoll, Paulettte Tercier,Ana Abad, Dave Goertz and Dave Butler for their comments and assistance throughout Inaddition, I would like to thank Ken Whittall for his assistance in analysing the data.This project ran considerably smoother due to advice and assistance provided by thetechnicians, in particular Yvonne Dumas and Ray Rodway.Funding for this project was provided by Shell Canada Limited, Esso Resources CanadaLimited, Petro-Canada Resources, and the Natural Sciences and Engineering Research Council ofCanada.xiiCHAPTER ONEIntroductionThe physical properties of porous rocks (i.e. elastic wave velocities, electrical resistivity,dielectric constant) are all governed, to varying extents, by the microgeometry of the pores andpore fluids. Empirical laws are commonly used to link physical properties to characteristics ofthese components; however, such laws are generally quite limited in application. In order todetermine more fundamental relationships between physical properties, the pore space, and fluiddistributions, analytical methods are required which can provide accurate characterization of thepore space and fluid distribution of rocks.The methods currently available for characterizing porous media include mercury injectionporosimetry, gas adsorption/desorption methods, electron and optical microscopy, and nuclearmagnetic resonance (NMR).Mercury injection porosimetry and gas adsorption/desorption methods operate on the samebasic principle: fluids which are not attracted to the surface of rocks must be forced into a rock byincreasing the pressure on the fluid. The size of the pores or pore throats through which the fluidenters or escapes from the rock is related to the fluid pressure in the system. The results of thesetwo methods are generally presented as plots of the volume fraction of the total void space which isaccessible through specific of pore throat sizes.The most accurate characterization of porous media is accomplished using microscopy toimage the pore networks. Both electrical and optical microscopy are used to image the surfaces ofthin or polished sections of rocks. Pore sizes may be determined from microscope images eitherthrough statistical extrapolation of three-dimensional pore sizes from two-dimensional images ordirect measurement of pore sizes from a three-dimensional image of the pore network. Highquality, two-dimensional images have been collected for many years; however, pore networks areinherently three-dimensional. Thus one of the goals of this thesis was to attempt to collect two-1dimensional images which could be used to reconstruct a three-dimensional image of a porenetwork.Three-dimensional data have been obtained through the collection of stereo pairs with thescanning electron microscope (SEM) and optical microscopes (Pitmann and Duschatko, 1970). Inaddition, three-dimensional images have been reconstructed from vertical series of two-dimensional images, or serial sections, produced through physical sectioning of the media (Straleyand Minnis, 1983, and Lin and Hamasaki, 1983). The method used in this thesis avoids many ofthe problems associated with stereo pair analysis and physical sectioning of rocks by collecting aseries of two-dimensional images of a pore cast through optical sectioning of the sample with theConfocal Laser Scanning Microscope (CLSM).The CLSM was used to collect independent two-dimensional images as well as verticalseries of two-dimensional images of three samples: Berea 100, a quarry sandstone, and twocarbonates, samples S55 and S77. The independent two-dimensional images were collected todetermine the pore sizes of the samples. The serial sections were collected to determine if theCLSM could be used to produce three-dimensional images of a pore network.Nuclear Magnetic Resonance was the second method used to characterize porous media inthis thesis. Advantages of NMR include speed, non-destructiveness, and the fact that samples maybe wet at the time of the measurement. NMR may be used to determine pore sizes because thespin-lattice or longitudinal (T1) relaxation times are dependent on the geometry of the pore space.The relation between pore sizes and T1 relaxation times is based on the increase in the rate ofrelaxation in fluids near a surface compared to the relaxation rate in bulk fluids. This results in thespin-lattic relaxation of fluids in porous media being a function of the surface to volume ratio of thepores (SfVre).Prior investigations of NMR have focussed primarily on determining the porosity,permeability, and pore sizes of sandstones. The determination of fluid distributions and theinvestigation of pores in carbonates, which are irregularly shaped and often very large, have rarelybeen considered. Thus, the other objective of this thesis was to test the current theory of NMR on2carbonates, and to determine what information on fluid distributions may be obtained using thistechnique.NMR data on twenty-two carbonate samples, nineteen of which were from the sameformation, was collected and analysed using two approaches: a fast diffusion method and a slowdiffusion method. The fast diffusion method is theoretically applicable only to those sampleswhich contain very small pores (<20 microns) while the slow diffusion method is applicable tosamples containing pores of any size. The nineteen samples from the Mississippian ElktonMember of the Turner Valley Formation had very little visible porosity, and therefore it wasexpected that these samples would be within the fast diffusion regime. However, the porosity ofthe three samples from the Mississippian Turner Valley Formation, S54, S55 and S77, wasdominated by visible porosity and thus the accurate analysis of the NMR data for these sampleswas expected to be possible only with the slow diffusion method.Determining fluid distributions and the pore sizes of carbonates from NMR data has rarelybeen attempted and the CLSM has never been used to image pore sizes in rocks. However, theadvantages associated with these methods may result in more accurate characterization of the porespace and fluid distributions in rocks than has previously been possible.3Characterization of Porous Rockswith the Confocal Laser ScanningMicroscopeINTRODUCTIONThe physical properties of porous rocks (eg. elastic wave velocities, electrical resistivity,dielectric constant) are all governed, to a varying extent, by the microgeomeiry of the pore space.Empirical laws are commonly used to link physical properties to characteristics of the pore space,but such laws are generally quite limited in application. In order to determine more fundamentalrelationships between physical properties and the pore space, analytical methods are needed thatcan provide accurate characterization of the pore space in rocks. The objective of this study was todevelop a means of obtaining and analyzing three-dimensional digitized images of the pore spaceusing a confocal laser scanning microscope.Characterizing Porous MediaThe current methods used to determine pore network characteristics in geologic samplesinclude mercury injection, gas adsorptionldesorption, optical microscopy, and electronmicroscopy. Each of these techniques, and their associated problems, is briefly described below.4Mercury InjectionThe mercury injection method (Wardlaw, 1976) produces a size distribution which is morerepresentative of the pore throats than of the pores. The general procedure for mercury injection isto remove all fluids from the sample, and then, with the sample under vacuum, force the mercuryinto the evacuated sample. The mercury is forced into the sample by progressively increasing thepressure on the system. The volume of mercury entering the pores and the corresponding pressureincreases are recorded and plotted as a pressure-volume or mercury injection capillary-pressurecurve. This curve is used to calculate pore throat size distributions based on the assumptions thatthe pressure required to force mercury into a pore is related to the size of the pore throat. Thus, thevolume of mercury which enters the sample between specific pressure limits is related to thatfraction of the total pore volume which is connected by pore throats within specific size limits.The pore throats are frequently assumed to be tubular and the resulting equation relatingpressure to the pore size is:D _4acos9IC—where Pc is the pressure, a is the surface tension of the liquid (mercury), e is the contact angle ofthe liquid on the solid, and d is the diameter of the tube, or pore throat (Wardilaw, 1976).Gas Adsorption/DesorptionThe gas adsorption/desorption method described by Pierce (1953) is similar to the mercuryinjection method in that the pressure on the system is altered to either push the gas into the sampleor allow the gas to escape. An adsorption or desorption isotherm of pressure versus volume of gasadsorbed or desorbed is plotted and one of the isotherms is used to calculate the pore sizes of the5sample. The isotherm selected to calculate the pore sizes is always the one which representscapillary equilibrium in the sample. This is checked by comparing the surface area which can becalculated from the isotherm with the surface area calculated using the BET method (Brunauer etal., 1938). For either isotherm, the volume at the end of the pressure step is assumed to be relatedto the size of the pore which is filling or emptying, and, using the Kelvin radius computation, theradius of the core within the pores which is not filled (rk) may be determined:r= 4.14log(j)where P0 is the initial pressure and p is the current pressure. The film of gas on the walls of thepores is assumed to contain “n” statistical layers and thus the pore radius (rn), is equal to:r = rk+3.6nwhere 3.6 A is the thickness of an adsorbed layer of nitrogen.Electrical and Optical MicrosconvThe most accurate characterization of pore networks in three-dimensions is based, as allpore characterization was in the early twentieth century, on the use of microscopy to image thenetworks. Both electrical and optical microscopy are used to image the surfaces of thin or polishedsections of rocks. Currently, most researchers are focussing on two objectives: producing highquality images which can easily be separated into the pore and grain components, and collectingand analyzing three-dimensional data.High quality, clear images are often produced by injecting a material into the pore spacewhich will enhance the contrast between pore space and grains in the rock. Conventionally, rock6samples in which the voids are to be examined were injected with a blue-dyed epoxy prior to thinsectioning. Other methods of improving the contrast between the pore and solid include injectingthe pore space with otherwise modified epoxy, woods metal, or wax (Straley and Minnis, 1983,Soeder, 1990, Winnsauer et al. 1952, and Nutting, 1929). These methods all make it possible tosee pores which, due to their small size or proximity to a large grain, would otherwise beimpossible to see. However, even using this technique, microporosity is often difficult todifferentiate from the surrounding grains because many of the injected materials are translucent inthe small quantities required to fill such porosity.During the past two decades, the imaging of microporosity has been achieved using boththe fluorescent optical microscope and the scanning electron microscope (SEM). Fluorescent dyeshave been employed to render microporosity visible to fluorescent optical microscopes (Soeder,1990, Ruzyla and Jezek, 1987, Rassineux et al., 1987, and Gies, 1987) and epoxy pore casts havebeen used to create high contrast images using the SEM (Pitmann and Duschatko, 1970, Slraleyand Minnis, 1983, Lin and Hamasaki, 1983). The pore casts examined with the SEM consist ofuntreated epoxy casts (Pitmann and Duschatko, 1970) and “double” pore casts, in which anuntreated and a brominated epoxy were used. These “double” pore casts were produced byinjecting epoxy into the pore space, dissolving the rock with hydrochloric and/or hydrofluoric acidand then injecting a second epoxy in place of the rock (Straley and Minnis, 1983, Lin andHamasaki, 1983). Due to the large difference in atomic number between untreated and brominatedepoxy, a high contrast back- scattered SEM image is obtained from such double pore casts,allowing very small pores to be imaged.The methods just described solve the problem associated with the low contrast between thepore space and grains in most geological samples; however, the lack of quantifiable threedimensional data is a problem which, until recently, was not addressed. Three dimensional data isimportant because pore networks are inherently three-dimensional while most images are only twodimensional. Three-dimensional data has been obtained both with the SEM and opticalmicroscopes by collecting stereo pairs of pore networks (Pitmann and Duschatko, 1970); however,7these images are very difficult to interpret quantitatively. The use of automated computer analysisis making this approach more viable, as Lin and Hamasaki (1983) have shown in their work onthree-dimensional reconstruction and analysis.The only other approach which has been attempted, in order to obtain three-dimensionalimages, was to physically section the sample to obtain a series of two-dimensional images, orserial sections. Two methods routinely used to produce serial sections are: successively grindinglayers off the surface of a sample (Lin and Hamasaki, 1983), and producing a two-epoxy porereplica which is then microtomed to produce successive thin sections of the epoxy (Straley andMinnis, 1983, Lin and Hamasaki, 1983). Both of these methods allow a “reconstruction” of theseries of two-dimensional sections to fomi a three-dimensional picture of the network. Thedrawbacks associated with these methods are primarily due to the distortion of the samples duringthe physical sectioning and the fact that large amounts of time and skill are required to produce thedata. When a sample is ground, the exact amount removed between imaging must be the sameeach time and, more importantly, the successive surfaces must be parallel, If distortion of theslices is to be minimized, care is also required when slices are produced with a microtome. Evenwhen a microtome is skillfully operated, a permanent distortion of five percent is usually expectedand the force of the blade used to cut the sample may cause the sample to twist, thus introducingmore error (Straley and Minnis, 1983).The method used in this thesis to image pore networks avoids many of the problemsencountered with the methods described above. A pore cast of the sample can be used to produce aseries of two-dimensional images without physically sectioning the sample, and the collection andanalysis of the images is relatively quick and easy. Also, the pore cast is not altered in any way,therefore re-examination of the same areas in their original orientation is always possible. Theproduction of serial sections without physical sectioning is made possible by the use of theconfocal laser scanning microscope (CLSM), the history and capabilities of which are describedbelow.8THE CONFOCAL LASER SCANNING MICROSCOPE (CLSM)The CLSM is one of many microscopes made possible by the introduction of theconventional scanning microscope (Young and Roberts, 1951). Conventional microscopes imageover a large field of view, while scanning microscopes are only required to image one point at atime and the final image is constructed through scanning. Removing the requirement to image alarge area allowed the introduction of modified optical systems for specialized use. One of thesemodifications produced the confocal scanning microscope.The commercial confocal scanning microscopes available now are generally reflectionand/or fluorescence microscopes. The light source is a laser focussed through a pinhole by acondenser lens. The light from this point source is reflected through the objective lens by adichroic mirror which focusses it on the sample. Reflected or fluorescent light from the sample iscollected by the same objective lens, passed back through the dichroic mirror and focussed into adetector pinhole. The detector, often a photomultiplier tube, measures the intensity of the lightreflected or emitted by the sample (Bacallao and Steizer, 1989).The confocal optical system described above was first described by Minsky (1957), whorecognized its possibilities for enhanced resolution and depth discrimination. The enhancedresolution is a result of Lukosz’s principle, which states that resolution may be increased at the costof the depth of field (Lukosz, 1966). The depth discrimination or optical sectioning property of theconfocal microscope is produced by the combined use of a point source and point detector. Thus,light from either above or below the focal plane is focussed behind or in front of the detectorpinhole and is not imaged. The accuracy of the depth discrimination has been discussed andcalculated theoretically by many researchers: for a dry lens with 0.95 numerical aperture, the depthof field is approximately 0.24 .tm and for oil immersion lenses, even smaller depths of field arepossible (Russ, 1990).The loss of a large depth of view makes the confocal microscope unsuitable for thoseapplications where a three-dimensional view is required or the sample is very uneven, causing a9great deal of scattering and thus degrading the image. However, there are many applications inwhich the optical sectioning property is actually a requirement for producing high quality images.Applications of microscopy which have been made possible by the optical sectioning of theconfocal microscope include many previously iniractable problems in the biological sciences(Wilson, 1990 and Shuman et al., 1989). The primaiy difficulty facing most biologists usingconventional microscopes was the haze of out-of-focus fluorescence produced by overlappingtissues. This haze hopelessly degraded the images and was first addressed by microtoming cellswhich had been fixed in plastic or wax. Slices 2 lim to 10 p.m thick were removed with themicrotome, transferred onto a slide, and examined with a microscope (Marshall et aL, 1990).Microtomoscopy, however, has similar difficulties in biology to those encountered in geology. Animmense number of slides have to be imaged, it is difficult to re-align the images, and distortion inthe original architecture of the cell are caused both by the fixation process and by microtoming.The majority of the difficulties encountered in mechanically sectioning cells are avoided bythe use of the CLSM. Cells of any shape and size may be viewed, and most importantly, they maybe viewed without fixation. In some cases, they may be imaged while functioning (Bacallao andStelzer, 1989, Egger and Petran, 1967). Thus a cell process, from neuron activity to bacterialgrowth, may be monitored as it occurs (Masters et al., 1991). These applications are only a smallportion of those possible with the confocal microscope. Recently, a CLSM has been used inmetallurgy to count pores in sintered alumina (Russ et al., 1989), in engineering to do surfaceprofile measurements (Hamilton and Wilson, 1982), and by physicists to image the surfaces ofmicrochips (Shuman et al., 1989). The success of the scanning confocal microscope in thesefields suggest that it could also be an effective technique for characterizing porous media in theearth sciences.A more complete treatment of possible applications and adaptations of the confocalmicroscope for specialty work is contained in Confocal Microscopy (Wilson,1990).10EXPERIMENTAL PROCEDURESSample PreparationTwo sedimentary rock samples, one high porosity, high permeability quartz sandstone, anda low porosity, low permeability carbonate, were used to develop and test the epoxy impregnationtechniques employed to create pore casts for this study. (Crocker et al., 1983). The quartzsandstone, Berea 100, has frequently been used in studies of porous media, therefore it wasselected to determine the accuracy of the imaging and image analysis techniques. The carbonate,referred to as samples S55, is Devonian in age and composed primarily of calcite with less thenone percent dolomite and silicate minerals. The low permeability and low porosity of sample S55increased the difficulty of thoroughly impregnating this sample with epoxy. Thus the efficiency ofthe impregnation techniques tested were judged by how thoroughly S55 was infiltrated.The samples were cut with a rock saw to expose two parallel surfaces and small particleswere cleaned off using compressed air. The dimensions of sample S55 was approximately 3mm x4mm x 4mm. The Berea 100 samples were disks 8.9 mm in diameter and 3 to 4 mm thick. Thecarbonate sample was smaller than the Berea 100 sample due to limitations in the available sample.A pore cast of the third sample considered in this study was prepared by Core Laboratoriesin Calgary; the others were not prepared by Core Laboratories because their techniques were notsuited to impregnating samples with permeabilities as low (< lmd) as that of sample S55. Thisthird sample is a high porosity, high permeability limestone referred to as sample S77. S77 isDevonian in age and composed primarily of calcite with approximately ten percent dolomite andless than one percent silicate minerals present.11MaterialsImpregnation of permeable media with epoxy and/or other fluids has been used in suchvaried fields as soil science, biology, and geology. Biologists use epoxy to fix tissues in place forexamination with microscopes. Soil scientists employ epoxy and wax to stabilize soils so the interrelationship of voids and grains in the samples may be studied. Geologists use epoxyimpregnation to increase the visibility of the void space in their samples; however, a geologist’srocks are sturdier than soils and larger than tissue samples, thus different epoxies and techniquesare used to impregnate these three types of samples.Both the characteristics of the sample and the purpose for which it will be impregnatedmust be considered when selecting an epoxy. Many rock samples considered have a very lowpermeability, thus a low viscosity epoxy is required. The formation of a pore cast requires theexposure of the epoxy to dilute hydrochloric and/or hydrofluoric acid, therefore it must not reactwith these fluids once it is cured. Finally, the low permeability of many rock samples requires thatthey be immersed in the epoxy for several hours to ensure complete impregnation. This means thatthe pot life, or the amount of time for which the epoxy retains its initial viscosity, must be relativelylong, with epoxies which require heating in order to cure being preferable to those which curerapidly at room temperature.There are at least two epoxies which match these requirements: Maraglas 655 resin with 555hardener (Soeder, 1990) and Spurr low viscosity embedding media (Spurr, 1969). Maraglass 655is a low viscosity epoxy resin with a pot life of several days. The Spurr low viscosity embeddingmedia was introduced by Spurr (1969) and has a viscosity of 7.8 cP at 20°C, similar to that ofwater, which is 1.002 cP at 20°C. Several different mixtures of the hardening and epoxycomponents are possible to produce fluids with different pot lives and curing schedules for theSpurr mixture. These possibilites are presented in Table 2.1. The mixture selected produced anepoxy with a pot life of approximately seven days and is the longer pot life/lower viscosity version12TABLE2.1SpurrLow-ViscosityEmbeddingMedia*StandardmediumSuggestedmodificationsofthemediumABCDEINGREDIENTFIRMHARDSOFTRAPIDCURELONGERPOTLIFELOWERVISCOSITYVinylcyclohexeneoxide(VCD)10.0g10.0g10.0g10.0g10.0gDigycidylether ofpolypropyleneglycol(DER736)6.0g4.0g7.0g6.0g6.0gNonenylsuccinicanhydride(NSA)26.0g26.0g26.0g26.0g26.0gDimethylaminoethanol(DMAE)0.4g0.4g0.4g1.0g0.2gCureSchedule(hr)at 70°C888316PotLife(days)3-43-43-427*Spuff,1969C)listed in Table 2.1. Both the Spurr and Maraglas epoxies could have been used to produce porecasts; however, Spurr’s epoxy resin was easily available and had a longer pot life, therefore it wasselected.The pore casts were to be imaged using the fluorescence mode of a CLSM, so the epoxyhad to be fluorescent before it was used to impregnate the rock samples. This requirement was notdifficult to meet as most epoxies, including Spurr low-viscosity embedding media, are naturallyfluorescent in the 460 to 490 nm range (Cercone and Pedone, 1987). This natural fluorescence isrelatively weak, therefore the fluorescence of the epoxy was enhanced by the addition of 0.1 g ofRhodamine B per 25 g of epoxy. The visibility of the pore cast to optical microscopes was alsoimproved by the addition of approximately 0.1 g of oracet blue dye per 25 g of the epoxy.Pore Cast ProductionThree impregnation techniques were tested to determine which method was best suited toour samples. In all techniques, the samples were dehydrated by heating in an oven at 105°C forseveral hours prior to impregnation. The dehydration was necessary because Spurr’s epoxy resinis incompatible with water, therefore the sample must be completely dry to enabe the epoxy to fullypenetrate the void space.The samples were impregnated and cured in a teflon sample holder. Teflon sample holderswere used because teflon is resistant to the 70°C temperature required for curing and it does notreact with either the epoxy or the two acids used to etch the samples. The holders were cylindricaland the sides tapered at approximately 5° so that a few taps on a counter would cause the sample tosimply slide out of the holder.After curing the sample at 70°C for at least 16 hours, the impregnated sample was removedfrom the holder and surface ground to expose the rock. The surface grinding was accomplishedwith grinding wheels, first with 250 mesh carborundum, then with 600 mesh (Al203)grit. A finalgrinding was done by hand with 900 mesh (A1203)grit on glass to smooth the surface. Problems14associated with grinding the sample include smearing the epoxy, removing too much of thesample, and not removing enough. Smearing of the epoxy over the grains and failing to removethe irregularities of the original surface may result in the pores appearing much larger than theiractual size. Alternatively, removing too much epoxy may result in the loss of the pore cast in caseswhere the epoxy did not fully penetrate the sample.The exposed surface was etched with acid to remove the rock, leaving a cast of the poresystem. Two acids were used: hydrochloric acid (HC1) to dissolve carbonate rocks, andhydrofluoric acid (F1F) to dissolve silicates.Impregnation TechniquesWet ImpregnationA variation of the ethanol replacement method used in soil science and biology (Marshall etaL, 1990 and Smart and Tovey, 1982) was attempted. In this technique, the dehydrated samplewas immersed in ethanol immediately after its removal from the oven. Several hours were allowedfor the ethanol to completely saturate the sample and the degree of saturation of the samples withethanol was checked by comparing the weight of the sample before and after soaking in theethanol. The volume of ethanol adsorbed was calculated and compared to the porosity of thesample; every sample was found to be at least 80% saturated with ethanol. The sample wasreturned to the ethanol for an hour before the process of replacing the ethanol with epoxy wasinitiated.Ethanol was replaced with epoxy by removing fluid (epoxy and/or ethanol) from thesample holder until it barely covered the sample, then epoxy was poured in to increase the fluidlevel to at least twice its original level. This was repeated every two to three hours throughout theday for two days, for a total of seven or eight changes in the fluid. The epoxy is soluble inethanol, thus if the two fluids mix thoroughly between each addition of epoxy, the amount ofethanol remaining in the fluid should be less than one percent. To ensure that the mixing was as15complete as possible, a special sample holder was constructed in which the sample was suspendedon a porous platform and a stirring rod below the platform agitated the fluid continuously to aid inthe mixing of the ethanol and epoxy. This method allowed the epoxy to easily penetrate even lowpermeability samples; however, the small amount of ethanol remaining in the epoxy after it was setreacted with the hydrofluoric acid, distorting or completely destroying the pore cast. Therefore thistechnique is not suitable for samples which will be etched with acid.Vacuum Impregnation: Immersion TechniqueThe simplest technique used required the dehydrated sample to be immersed in the epoxyimmediately after removal from the oven. The sample, immersed in epoxy, was placed within avacuum dessicator and a vacuum applied to remove the air from the voids and allow the epoxy topenetrate them. The sample remained within the vacuum dessicator for four to five hours, the firsttwo to three hours under vacuum and the final two hours equilibrating to the surrounding airpressure. Unfortunately, a high pressure vacuum was not available to pull the air out of the pores,thus the low permeability samples were not impregnated beyond the surface layer of the sample.Vacuum Impregnation: Lowering TechniqueThe most successful technique used a form of vacuum impregnation that required thesample and the epoxy to be degassed separately, then the epoxy either poured over the sample orthe sample lowered into the epoxy. The epoxy may be poured over the sample either from a beakeror a separatory funnel. The beaker must be suspended over the sample and it must be possible tomanipulate it without opening the vacuum chamber. Vacuum apparatus similar to this iscommercially available from several companies; however, a simpler solution is to use a loweringapparatus. This requires a vacuum vessel sealed with a rubber stopper through which a rod maybe extended into the vessel. The sample holder containing epoxy is placed on the base of the16vessel, and the sample is suspended over the sample holder by attaching it to the rod with string ornetting. The vacuum is applied to both the sample and the epoxy for an hour to remove all the airfrom both. The sample is then lowered into the epoxy and the netting or thread is released. Thevacuum is continued for approximately thirty minutes to remove any gas which is trapped in thevoids of the sample.EtchingFollowing the curing of the epoxy, the sample was removed from the sample holder andsurface ground to expose the rock surface as described previously. The exposed rock and epoxysample was immersed in five percent hydrochloric acid (HC1) for between five and sixty minutesdepending on the depth to which etching was desired. This dissolved any carbonate rock presentin the upper portion of the sample. Residual HC1 was removed by rinsing the sample repeatedlywith deionized water and then any silicate minerals present in the sample were dissolved byimmersing the sample in 5% hydrofluoric acid (HF) for between five and sixty minutes. The HC1was used first because HF reacts with carbonate rock to form insoluble calcium fluoride (Grayson,1956). Both HC1 and HF were tested to ensure they did not react with Spurr’s epoxy resin byPittman and Duschatko (1970). Residual acids were removed by repeated rinsing with deionizedwater before any imaging was done.All of the impregnation techniques described above worked quite well for the Berea 100sample; however, the microstructure of the pore network in the low permeability carbonate, S55,was too fragile and collapsed when the surrounding mineral was removed. This problem wasovercome by producing “double” pore casts of this sample.17Double Pore CastsDouble pore casts were produced by re-impregnating the sample with a different epoxyafter the carbonate portion was dissolved with HC1. The purpose of this second epoxy was tosupport the pore cast. The second epoxy used was London Resin White (LRW), a transparent,low viscosity, non-fluorescing resin with a very long pot life. This epoxy was cured at 72°C for atleast 24 hours. Impregnation of the sample with a second epoxy stabilized the microstructure ofthe void space enough to enable the subsequent removal of silicate minerals in the sample using 5%HF for between five and sixty minutes. The LRW was tested to ensure it did not react with eitherHC1 or HF by immersing a block of it in both fluids for several hours and subsequently examiningthe block microscopically. The examination revealed no alteration in the character of the LRW so itwas assumed safe to use in the production of these “double” pore casts.DATA COLLECTION AND ANALYSISImages were collected using an MRC BioRad 500, a highly versatile microscope which iscommercially available and operates either in reflection or fluorescence mode. The imagescollected with the BioRad 500 were 768 by 512 pixels in size and each pixel had an intensity valuewhich corresponded to the intensity of the fluorescence at that point in the sample, with 0 being noresponse and 255 being the most intense response. The portion of the sample being imaged wasscanned 10 to 15 times and the images were stacked to improve the signal to noise ratio. This wasaccomplished interactively; when the clarity of the image ceased to improve, the scanning washalted. The images were stored in a binary format, thus the algorithm written in the Cprogramming language and listed in Appendix 1 was required to transform them into a format(ASCII) readable by the FORTRAN programs written to analyse the images.The primary requirement for the collection of data is that enough images be collected toprovide a statistical sampling of all pore and pore throat sizes. Forty two-dimensional images of18Berea 100, twenty-one of sample S55, and thirty-three of sample S77 were collected. Some ofthese images were part of a vertical series of contiguous images, and therefore only one of theseries was analysed to determine pore sizes to avoid biasing the pore size distribution towards thepores in that portion of the samples. The final analyses were accomplished with sixteen Berea 100images, twenty-one S55 images, and thirteen S77 images. These images were manipulated toproduce pore and pore throat distributions using the process described below, which wasimplemented using the programs listed in Appendix 2.Initially, the image was altered from one with an intensity range of 256 to one which wasbinary. The process of altering an image with a large intensity range to a binary image is oftenaccomplished by visual means using an image analysis program or by using the intensity contrastbetween adjacent pixels to choose the boundaries of the objects. The separation of these imageswas accomplished similarly by using an imaging program to select the intensity ranges whichcorresponded to epoxy and void (pores and grains). Also, several rows of each image wereplotted and the separation value chosen to separate pore from solid was super-imposed on this plot.A value which lay below the peaks which represented the pores and above the lows whichrepresented the solid regions was chosen. Examples of the value chosen super-imposed on a plotof the image intensity are shown in Figure 2.1 for two rows in the image referred to as ABR1.The resulting binary images when the separation value chosen is varied are shown in Figure 2.2.The small holes in the pores and solid which may have been due to either the presence of verysmall pores or to uneven fluorescence intensity in the epoxy were removed. This wasaccomplished using a dilation/erosion sequence followed by an erosion/dilation sequence. Theseare standard procedures in image analysis. To illustrate the procedure, the process of removingsmall pores within the solid is described. First the solid is dilated by altering any pore pixels incontact with a solid to a pseudo-solid, then any pixels added to the solid via dilation which were incontact with a pore are changed back into pore pixels. Thus, a pore which was less than 1.5 pixelsin diameter or, for the images collected at the highest magnification, 3.0 tm in diameter, waschanged to solid. Similarly, small solid areas completely within a pore were removed by first19IMAGE ABR1 - ROW 40>zzPIXEL NUMBERIMAGE ABR1 - ROW 250I I I30020010003000 200 400 600 800200l000•LLLrJ - fiJFJ J IrU\%J%0 200 400 600 800PIXEL NUMBERFigure 2.1: Pixel intensities in rows 40 and 250 of image ABR 1. The intensity value 50 isindicated as possible separation value between pore pixels and solid pixels.20Figure 2.2a: Binary versions of Berea 100 image ABR1 with pore/solid separation valueof 40. Black = solid, and white = pore space.2OOn212OO.zmFigure 2.2b: Binary versions of Berea 100 image ABR1 with pore/solid separation valueof 50. Black = solid, and white = pore space.222OOiniFigure 2.2c: Binary versions of Berea 100 image ABR1 with pore/solid separation valueof 60. Black = solid, and white = pore space.23dilating the pores, then removing (eroding) those pixels which remained in contact with the solid.The pores were separated and those points within the pores which were farthest from theedges of the pores were found using the watershed algorithm described by Russ (1990). The mostcentral coordinates of each pore and pore throat were recorded so that the size of the pores and porethroats could be determined.The size of the pores and pore throats was determined using methods similar to thoseemployed by Ehrlich et al. (1984) and Doyen (1988). The pixels within the pores and pore throatswhich are most central are used to compare the pore or pore throat to a previously constructedcircle in order to determine the largest diameter circle which is completely enclosed by the pore orpore throat.The radii measured in the manner described above are presented as volume fractiondistributions which are the ratio of the area of the circles described by each individual radius to thetotal area of all of the measured pores. This ratio is assumed to be proportional to thecorresponding three-dimensional ratio of the volume of a pore to the total volume of the pore spacebecause, if sufficient sections are randomly imaged from a sample, the sizes, numbers, and areafractions of the pores images should be equal to the sizes, numbers, and volume fractions of thethree dimensional objects (Wicksell, 1926, and Russ, 1990).RESULTSThe images produced with the confocal laser scanning microscope were clear, high-contrastimages which clearly show a major variation in pore shapes and sizes between the three samples.The Berea 100 images (Figure 2.3) contained the expected network of larger pores connected bymuch smaller pore throats. The carbonate images were dominated by thin, plate-like pores whichwere often associated with pores greater than 100 urn in diameter (Figure 2.4a). S55 containedonly a few of the larger pores (Figure 2.4b); however, S77 contained a large number of such pores(Figure 2.5).242OOimFigure 2.3: Original version of a confocal laser scanning microscope image of Berea 100.The image is called ABR1. The pore space is the brightest areas in the image.252OOimFigure 2.4a: CLSM image of vugs in the carbonate sample, S55. The image is called F78and the brightest areas of the image represent porosity.26Figure 2.4b:CLSM image of vugs in the carbonate sample, S55. The image is called T78and the brightest areas of the image represent porosity.200qxn27Figure 2.5: Original version of a CLSM image of carbonate sample S77. The image is called Fl.The pore space is the bright areas in the image.2829The choice of a “cut-off’ or pore/solid separation value to alter the original images to binaryones was different for each sample. The magnitude and overall distribution of porosity changedvery little for the Berea 100 image, referred to as ABR1, for solid/pore separation values between40 and 60 (Figure 2.2), therefore the final value selected was 50. This value was tested on theother Berea images and found to be appropriate for them as well.The S55 images were more difficult to separate into two phases. The fluorescence of thethin, plate-like pores was faint compared to the larger pores, thus when a high separation valuewas chosen, the thin pores were neglected (Figure 2.6a), while if a lower value was used, theresulting image contained porosity which, in the original, seemed to be an out-of focus haze(Figure 2.6b).Separation of the S77 images into pores and solid was relatively easy as the intensity of thepore and solid regions were well separated in these images (Figure 2.7). The pores and porethroats were measured from the coordinates which were farthest from the edges of the pores orpore throats (Figure 2.8). Pore size and pore throat distributions were obtained from the binaryimages of Berea and S77, and from those binary S55 images which were produced using a highseparation value. The pore size distributions presented here are volume fraction distributions,therefore the percentage associated with each pore radius indicates the fraction of the total porevolume occupied by pores of that size.The Berea 100 pore radius volume fraction distributions all have a similar shape (Appendix3), as shown by the examples in Figure 2.9a. The distributions extend from 1 to 60 $.tm with themaxima occurring primarily between 10 and 30 im with several isolated peaks above 30 I.Lm.Thepore size distributions of several of the images combined also extends from 1 to 60 .tm and themaximum occurs between 10 and 20 Lm (Figure 2.9b).The pore size distributions for sample S55 were quite different (Figure 2.lOa) and a copyof the distribution for each image is contained in Appendix 4. These distributions were lessuniform, with most maxima occurring between 7 and 20 Lm and a few between 20 and 60 I.tm.The combined pore size distributions of all of the images extends between land 50 Im and reaches30Figure 2.6a: Binary version of carbonate S55 CLSM image F78 (Figure 2.4a) produced using asolid/pore separation value of 50. White = porosity, and Black = solid.2OOun31Figure 2.6b: Binary version of carbonate S55 CLSM image F78 (Figure 2.4a) produced using asolid/pore separation value of 100. White = porosity, and Black = solid.200pmw ‘32Figure 2.7: Binary version of carbonate S77 CLSM image Fl (Figure 2.5) produced using asolid/pore separation value of 15. White = porosity, and Black = solid.33C)IFigure 2.8: Binary versions of Berea 100 image ABR1 with pore throats separated from thepores. The centers from which the sizes of the pores and pore throats will be measured areindicated by a contrasting colour within the void space. Black = solid, the White = porosity, andthe dotted areas = pore throats.35Figure 2.9a: Examples of the Berea 100 pore size distributions. The upper plotrepresents the pore size distribution of image ABR1, the lower of image HBR1.Image ABR11060.—I.J0C)010Pore Radius (m)0.120.080.040.000.080.060.04 -0.02 -0.0010Image HBR1&I ii106 10 -IPore Radius (m)370.—00.085 ImagesA0C.)J0.060.040.0210Pore Radius (m)Li0416 Images0.080.060.040.020.00Figure 2.9b: Berea 100 cumulative pore size distributions. The upper plotrepresents the summation of the distributions of 5 images, the lower of 16 images.Pore Radius (m)380.-40000Image F780.30.20.1•0.01060.20.110Pore Radius (m)i04Image N780.010L10Pore Radius (m)Figure 2. 10a Examples of the carbonate (S55) pore size distributions. The upperplot represents the pore size distribution of image F78, the lower of image N78.3920 ImagesFigure 2. lOb:Carbonate (S55) cumulative pore size distributions. The upper plotrepresents the summation of the distributions of 5 images, the lower of 20 images.5 Images0C)J00Pore Radius (m)i040.120.080.040.00 -i106ii10Pore Radius (m)i0440a maximum at approximately 7 pm (Figure 2. lOb).The S77 pore size distributions (Appendix 5) resembled discrete pore size spectra , with thelargest pores present in the image dominating the distribution (Figure 2.11 a). However, the largestpores present in the image varied from 50 to 220 pm in radius, resulting in a cumulativedistribution which extended primarily between 10 and 220 pm with several maxima between 70and 150 p.m (Figure 2.1 lb).The mean values of the volume fraction pore size distributions were calculated and arecollated in Tables 2.2, 2.3 and 2.4 for the Berea 100, S55 and S77 samples respectively. Themean pore size of the Berea images varied from 20.2 to 30.4 p.m, with an overall mean pore size of25.3 microns and a standard deviation of 15 pm for all of the images considered together. Thepore sizes in S55 are generally smaller, ranging from 8.9 to 33.0 pm, with an overall mean poresize of 17.5 p.m with a standard deviation of 11 pm. The largest pore sizes and largest variation inmean pore sizes occur in sample S77. The mean pore size of this sample ranges from 30.6 to150.9 pm for the individual images with an overall mean pore size of 90.2 p.m and standarddeviation of 54 p.m for the cumulative pore size distribution.The porosity of the Berea 100 images varied a great deal, from 0.26 to 0.41 (Figure 2.12),with the majority of the images having porosities between 0.29 and 0.34. There is a similarvariation in the porosity of the images for the carbonates, S55 (Figure 2.13) and S77 (Figure2. 14).The minimum and maximum porosities of the S55 images produced using a high separationvalue are 0.04 to 0.16, with a single outlier at 0.25. The porosity of these same images when alow separation value was selected was much higher, ranging from 0.11 to 0.39, with most of theimages having porosites between 0.15 and 0.25. S77 image porosities varied from 0.02 to 0.18with an average of 0.088.410.—a)E00.-4C)0Image JiImage Gi4AA1l106 io4Pore Radius (m)0.50.40.3-0.10.00.4_0.30.1 -0.0- I106 io’ io3Pore Radius (m)Figure 2.1 la: Examples of the carbonate (S77) pore size distributions. The upperplot represents the pore size distribution of image G 1, the lower of image Ji.I425 Imagesio-5Pore Radius (m)13 Images0.20.10.010.080.060.040.020QC>0C-)0) ‘.3liii I0.001O..6 10”Pore Radius (m)Figure 2.1 ib: Carbonate (S77) cumulative pore size distributions. The upper plotrepresents the summation of the distributions of 5 images, the lower of 13 images.43C’)LULL0LUDzBEREA IMAGES - POROSITYPOROSITYFigure 2.12: Porosity (the ratio of pore pixels to the total number of pixels) of theBerea 100 images.CON. CD O CD.- CJ CD L() CON. CD O —0000 000000000 044S55 IMAGES - POROSITY2_____LOW SEPARA11ON VALUEMEAN POROSITY: 0.220 - . —v—v—. . . - . . - . .i— .. •— -C’ e Ifl %O N 00 O- ‘f ‘D N 00 O\ -C4 fl ‘.0 N 000’.POROSITYS55 IMAGES - POROSITY3 HK3H SEPARATION VALUEMEAN POROSITY: 0.11C,,10Il’. ‘.0 N 00 0’ — - (‘. ‘.0 N 00 0’. - C.l e - In00 0 0 0 0 0 0 0 0 0 0 0 0 0 0POROSITYFigure 2.13: Porosity (the ratio of pore pixels to the total number of pixels) of thecarbonate S55 images. The upper distribution represents the binary imagesproduced using a low intensity separation value, and the lower distributionrepresents those images produced with a high intensity separation value.45NUMBEROFIMAGES0•1%)0.020.03___________________Ci)0.040.05________—0.06Z0.070.08C-)0.09I.<0o.iCI)CD‘—0.11_0.12C-0.13o_________0.14CD0.15Ci)—0.160.170.18___________________________________CD 00)TABLE 2.2 Berea 100 Mean Pore SizesImage Mean Pore Standard Number of Mean Pore StandardRadius Deviation Images Radius Deviationjim jim jim jimABR1 24.8 13 1 24.8 13BBR1 23.7 15 2 24.2 14CBR1 22.8 13 3 23.7 14DBR1 26.4 16 4 24.4 15FBR1 20.2 9 5 23.7 14GBR1 21.5 12 6 23.3 14HBR1 26.4 17 7 23.8 14JBR1 28.6 19 8 24.3 15KBR1 26.1 17 9 24.5 15LBR1 24.9 13 10 24.6 15MBR1 30.4 17 11 25.2 15NBR1 27.4 17 12 25.4 15OBR1 22.0 12 13 25.2 15PBR1 24.3 14 14 25.1 15QBR1 22.6 10 15 24.9 15RBR1 29.7 17 16 25.3 15TABLE 2.3 Carbonate (S55) Mean Pore SizesImage Mean Pore StandardRadius Deviationjim j.tmA78 12.2 6B78 18.4 11C78 17.2 10D78 15.7 10E78 15.7 10F78 10.6 6G78 14.4 8H78 8.9 4J78 19.0 10K78 14.0 6L78 12.2 6M78 14.1 8N78 17.1 13078 12.0 7P78 15.7 9Q78 29.2 16R78 12.1 7S78 18.0 9T78 21.0 9U78 33.0 13Number of Mean Pore StandardImages Radius Deviationj.tm jim1 12.2 62 17.0 113 17.1 104 16.6 105 16.4 106 15.8 107 15.6 108 15.0 109 15.5 1010 15.4 1011 15.2 912 15.1 913 15.4 1014 15.1 1015 15.2 1016 16.6 1117 16.3 1118 16.4 1119 17.0 1120 17.5 1147TABLE 2.4 Carbonate (S77) Mean Pore SizesImage Mean Pore Standard Number of Mean Pore StandardRadius Deviation Images Radius Deviationj.im jim jim jimFl 113.6 50 1 113.6 50Gi 150.9 71 2 126.3 60Hi 38.2 13 3 121.8 62Ii 30.6 13 4 119.3 63Ji 97.6 33 5 116.3 60Ki 50.5 14 6 111.3 60Li 99.2 48 7 108.1 58Ml 57.5 25 8 103.3 57Ni 35.7 19 9 100.7 5801 75.1 53 10 97.3 58P1 63.4 26 ii 93.6 56Qi 90.2 43 12 93.2 55Ri 56.1 26 13 90.2 54DISCUSSION AND CONCLUSIONSThe confocal laser scanning microscope produced very clear, easily analysed images whenthe porosity being imaged was relatively uniform in size. Unfortunately, the S55 carbonate sampleconsidered in this project contained pores of widely varying sizes, and the enhanced fluorescenceof the larger pores overshadowed that of the small, thin pores. Thus the S55 images, althoughqualitatively useful, were very difficult to convert to a binary format. The two-phase imagesproduced by using widely separated cut-off values (Figures 2.6a,b), when compared to theoriginal image (Figure 2.4a) illusirate the primary difficulty in analyzing these images: the thinpores are either neglected completely (Figure 2.6b) or they are represented as being much thickerthan they actually are (Figure 2.6a). Sample S77, although it also contained pores of widelyvarying sizes, did not produce the same problems. This may be a result of two effects: theepoxy/dye mixture used by Core Laboratories to impregnate this sample produced a strongerflourescence, or, considering Core Laboratories lack of success in impregnating S55, the smallestpores may not have been penetrated using their technique.48A possible solution to imaging samples with two or more widely separated pore sizes is toimage the samples at two magnifications. The higher magnification images may be used todetermine the size of the smaller pores and the low magnification images to correlate the number orvolume of small pores associated with the larger ones. In this way the final pore size distributioncould be determined by estimating the number and size of the associated smaller pores when thesize distribution of the larger pores is determined.The choice of which pore/solid separation value to use in further analysis of the S55carbonate images was based on the porosity of the resulting images. The measured porosity of thesample was 0.048 (Chapter 3), therefore the largest value of the cut-off reproduced the correctporosity most accurately with an average porosity of 0.11, (Figure 2.10), and these images wereused for all further analysis.The large porosities of the images, an average of 0.322 compared to the measured value of0.21 for Berea 100 (Tercier, 1992) and an average of 0.088 compared to the measured value of0.09 for S77 (Chapter 3), are primarily due to experimental bias. The images are collected byscanning the surface of the sample and “randomly” selecting a portion to image. The selectionwas accomplished interactively, thus areas with little or no porosity were not imaged. If thisexperiment were repeated, one possible improvement is to mechanize the imaging to ensure that theentire surface of the sample is imaged regardless of the porosity included in the image.Alternatively, the number of image areas containing no pores could be recorded and included in thecalculation of the total porosity of the sample.An example of the locations of the centers used to measure the sizes of the pores and porethroats in the images is present in Figure 2.8. This figure shows how these centers were spaced toensure a minimum amount of overlap between the areas included within the measured circles. Thisresults in elongated pores containing three or more centers, thus illustrating the primary difficultyin determining pore sizes: their shapes are not easily described. The assumption of circular poreswhich was used here is frequently used, despite its inaccuracy, because it simplifies mostcalculations. A more representative shape for pores would be elliptical; however, measuring the49sizes would be complicated by the irregularity in the outline of the pores. A possible solution tothis difficulty is to approximate the size of ellipse which may be contained within a single porebased on the circles already used to measure the pores. This could be accomplished bydetermining the length of the major and minor axes of the ellipse by summing the diameters ofcircles paralleling the major axis and averaging those paralleling the minor axis.The overall mean pore sizes of the Berea 100 and carbonate images are listed below (Table2.5). The results for the Berea 100 images were very encouraging, showing good agreement withmercury injection pore size distribution of Berea 300 presented by Crocker et al. (1983). TheBerea 300 pore size distribution may be compared to that of Berea 100 as they are both Bereasandstones, the only variation between the the two being that Berea 100 has a permeability ofapproximately 100 mD, while Berea 300 has a permeability of approximately 300 mD. The Berea300 distribution (Figure 2.15) has a mean pore radius of 9.0 jim, compared to a cumulative meanvoid radius of 25.3 p.m determined for the Berea 100 images. The 9.0 p.m determined byCrocker et al. (1983) is substantially lower, as expected since capillary presssure mercury injectiondata reflects the size of pore throats, making the mean size lower than the actual average pore size.TABLE 2.5 Mean Pore SizesSample Mean Void Radius Standard Deviationjim jimBerea 100 25.3 15S55 17.5 11S77 90.2 5450I40 30 20 10 0PoreDiameter(pm)Figure2.15:Berea300capillarypressureporesizedistributions.(fromCrocker etaL,1983)100BereaMean—18.0gmMode=10.0urn0.20.4124102040(31The primary objective of this study was to develop a method of collection and analysis ofthree-dimensional images of pore networks. Despite the lack of complete three dimensional imagesin this study, the confocal laser scanning microscope is a powerful tool which could be used toaccomplish this goal. Unfortunately, the fragility of the pore casts ensured that the series of two-dimensional images which were obtained (Figure 2.16) did not encompass a thickness greater thanthe average pore diameter and thus could not be used to directly determine three-dimensional poresizes. However, the correspondence between the Berea 100 pore size distribution and thatdetermined by Crocker et al. (1983) for Berea 300 indicate that CLSM imaging is definitely aviable technique for obtaining accurate images of pore casts.52C)’C)Figure2.16:Consecutiveopticalsectionsof Berea100.Thespacingwas9.6micronsbetweenthesections.Fromrighttoleft,theimagesarethesecond,thirdandfinal(fourth)imagescollected.ABR1wasthefirstofthisset(Figure2.3).Thebrightareasintheimagesarethepores,andthedarkerareasrepresent thesolid.CHAPTER THREECharacterizing Porous Rockswith Nuclear MagneticResonanceINTRODUCTIONThe characterization of porous media using Nuclear Magnetic Resonance (NMR) hasbecome a focus of research as scientists seek to use it to address various questions in manydifferent fields of endeavor. The initial application of NMR to describe the characteristics ofporous media occurred in the oil industry. The Nuclear Magnetic Logging tool (NML) was used tocorrelate the NMR response of petroleum reservoirs with the physical properties of the rock (e.g.porosity, permeability, free fluid index, irreducible water saturation). Recently, the more indirectapproach of determining pore sizes from NMR data and then calculating the physical propertiesfrom the pore sizes has been used. NMR has also been applied to the determination of sizedistributions in such diverse media as oil/water emulsions (Parker and Rees, 1972), biological cells(Brownstein and Tarr, 1979), and cheese (Callaghan et al., 1983).The objective of this study is to extend the uses of NMR in geology to the characterizationof carbonate reservoir rocks. The suitability of the NMR method for determining pore sizedistributions of carbonate rocks will be investigated using the “two-fractions in fast exchange”model and the more general theory of Brownstein and Tarr (1979). In addition, the potential useof NMR to determine the fluid distributions in carbonates will be examined.The significance of this research derives from the fact that the microgeometry of the porespace and the pore fluids controls the petrophysical properties (i.e. permeability, electricalconductivity, dielectric constant) of porous media. These properties are used in the petroleumindustry to interpret well-logging data and to estimate the capacity and production potential of oilreservoirs. The environmental industry uses the physical properties of rocks to predict the positionand movement of groundwater and contaminants through the groundwater system.54HISTORY OF NMR IN PETROLEUM GEOLOGYNuclear Magnetic Logging (NML)The basic principles of NMR were first published by Bloch (1946). The possibleapplications of this technique to the determination of reservoir properties were quickly recognizedby the oil industry, and the first patent for a well logging instrument based on the principles ofNMR was issued in the 1950s to the California Company (Jackson, 1984). The first nuclearmagnetic logging tools (NML) were used in the 1960s. The industry expectations of determiningfluid types, permeabilities, and the producible fraction of oil in place, were not met with thesetools. However, more recent NMR logging tools do provide these data (Neuman and Brown,1982). The unreliability of most of the NMR logging tools was due to difficulties with themagnetic field. Producing a strong, homogeneous magnetic field in the area suffounding theborehole where the electro-magnet must fit within the fifteen to twenty centimeter diameter of awell-logging tool was one of the problems. This problem has been addressed by the newtechnique invented in 1978 at the Los Alamos National Laboratory (Jackson, 1984). Thistechnique is based on the creation of a toroidal region of homogeneous radial magnetic fieldsurrounding the welibore at a specific distance from the center.The most recent development in the the field of nuclear magnetic well logging is thepresentation of a prototype magnetic resonance imaging log, MRIL (Coates et al., 1991). The datawhich can be determined with this tool include a “sand-size” porosity, a bulk volume ofirreducible porosity and a link to the permeability of the rock.” (Coates et al., 1991)55Relationships Between NMR and Physical PropertiesEmpirical relations between the physical properties of rocks and NMR data have beensuggested since 1956 when Brown and Fatt (1956) related the fractional wettability of oiffleldrocks to the free induction decay (FID) curve measured with an NMR apparatus. Correlations withthe porosity, permeability, irreducible water saturation, the movable fluid in the reservoir and theoil/water saturations were also found (Timur, 1969, Robinson, 1974, and Vinegar et al., 1989).Pore SizesThe indirect approach of using NMR to determine pore sizes and then estimating thephysical properties of the rocks was first proposed by Loren and Robinson (1970). The majorityof research on this technique for geological materials has taken place in the last decade using spin-lattice (T1) laboratory measurements of sandstones and calculating pore sizes using the “two-fractions in fast exchange” model (Brownstein and Tarr, 1977, 1979).The calculation of pore sizes from Ti data using the “two-fractions in fast exchange” modelis based on the assumption that the fluid in the pores is separated into two phases: a surface phasewith a short relaxation time and a bulk phase with a substantially longer relaxation time. Themajority of the nuclei are assumed to relax at the surface of the solid, an assumption which isaccurate only if the pores are small enough so that all nuclei in the fluid can diffuse to the surface ofthe pore in a time less than the average relaxation time of the pore. If these assumptions arecorrect, then there is a one-to-one correlation between pore sizes and relaxation times.Various techniques have been used to calculate average relaxation times or distributions ofrelaxation times which accurately reproduce the collected NN’ER data. These methods includematching the data using a set number of relaxation times, using a “stretched exponential”representation (Borgia et a!., 1990 and Kenyon et al., 1988), and using inversion methods toproduce discrete (Schmidt, 1986, Munn and Smith, 1987, and Gallegos et al., 1987, Whittall and56MacKay, 1989) or Continuous spectra of relaxation times (Gallegos and Smith, 1988, Gallegos etal., 1988, and Kenyon et al., 1989, Whittall and MacKay, 1989).The most recent research of Davies and co-workers (1990) has applied the theory ofBrownstein and Tarr (1979) in which the fast diffusion assumption is not required. This studyfocussed primarily on sandstones, as did most of the others mentioned above. This predeliction tothe study of sandstones is present because the pore network in sandstones is relatively uniform, themajority of oil reservoirs in the Unites States are sandstones, and the surface effect is morepronounced in sandstones than in carbonates (Loren and Robinson, 1970 and Schmidt et al.,1986).Fluid DistributionsThe microscopic distribution of fluid in porous media has a pronounced effect on theirphysical properties. This effect has been observed in laboratory measurements of compressionalwave velocities (Domenico, 1976, Domenico, 1977, and Knight and Nolen-Hoeksema, 1990,seismic attenuation (Bourbie and Zinzner, 1984), electrical resistivity (Longeron et al., 1989 andKnight, 1991), and dielectric constant (Knight and Nur, 1987). Although NMR is a logicalmethod for determining fluid distributions in sedimentary rocks, only two studies of partiallysaturated rocks have been completed with this goal in mind (Schmidt et al., 1986, and Davies et al,1990).THEORYAn excellent review of the basic theory of Nuclear Magnetic Resonance may be found inExperimental Pulse NMR. A Nuts and Bolts Approach (Fukushima and Roeder, 1981). Theportions of this theory which are relevant to this study are summarized below.57Various nuclei (i.e: ‘H, ‘3C) possess an intrinsic magnetic field, or magnetic moment (M).This is created by the interaction between the angular momentum, or spin, of the nucleus and theelectric charge associated with the atom.Let us consider nuclei which do not interact with each other, and thus represent isolatedmagnetic moments. These magnetic moments will be randomly oriented; however, if an externalmagnetic field (a3) is imposed, they will equilibrate to the field by assuming two possible lowest-energy positions: precession around the directions parallel and anti-parallel to the external field.The parallel orientation is favoured slightly by thennodynaniics, thus a cumulative magneticmoment (M0)parallel to the external field is produced. The precession frequency, also referred toas the Larmor frequency, is a function of the type of nucleus and the strength of the magnetic field.M0 is much weaker than the constant field, H0, therefore the change in the overall fieldcannot be accurately measured. However, if a radiofrequency (if) pulse operating at the Larmorfrequency is applied perpendicular to H0, the nuclei will absorb energy from this field and rotateaway from H0 and the rf field. The angle through which M0 rotates depends on the amount of timefor which the rf pulse is applied. Thus, a 90° or 180° pulse is a radiofrequency pulse whichcauses M0 to rotate through 90° or 180° respectively. Immediately following a 90° pulse, M0 isperpendicular to both the rf field and the external field and will produce an emf field whichsubsequently induces a measurable voltage in the if coil. Through this process of first producing acoherent magnetic moment, M0, by imposing an external field, and then shifting M0 90° from theexternal field by applying a radiofrequency pulse, the magnitude of M0 can be measured.Nuclei in condensed matter, solids or liquids, are coupled to the degrees of freedom,referred to as the lattice, of that medium. This spin-lattice coupling tends to cause the nuclearmoments and the lattice to evolve towards thermal equilibrium concurrently. Thus when theexternal field, H0, is applied, the evolution towards M0 requires that both the nuclei and the latticereach thermal equilibrium. The progression towards M0 from the original components, M, Mand M, is assumed to be exponential with two time constants: T, for the longitudinal componentM, and T2 for the transverse components, M and M (Figure 3.1).58MMFigure 3.1: Nuclear Magnetic Resonance Reference FrameH0M59Longitudinal Relaxation Time (T1)The longitudinal relaxation time is the exponential time constant representing the evolutionof the longitudinal magnetic moment (Mi) of the nuclei into the moment M3 when an externalmagnetic field, H0, is applied to the media. The equation representing this progression is:-Mz=*(Mz-Mo) [1]The evolution towards M0 from M is accompanied by an exchange of energy between thelattice and the spinning nuclei and so T1 is often referred to as the spin-lattice relaxation time. Themagnitude of the spin-lattice relaxation time T1 reflects the efficiency of the energy transferbetween the lattice and the spinning nuclei, thus T1 is a function of the type of atoms present andtheir interaction with the lattice.The spin-lattice relaxation time is usually measured using the inversion-recovery pulsesequence which consists of a 180° pulse followed by a 90° pulse. The 1800 pulse flips themagnetization around so that it is equal and opposite to the cumulative magnetic moment, M0. The90° pulse is applied a time (t) after the initial pulse and shifts that component of the magnetic fieldwhich is in the z direction to the axis perpendicular to the applied and rf fields so that its magnitudecan be determined. This is followed by several seconds in which no pulse is applied and then thetwo pulses are applied again with the time separation varying with each application until the fulldecay curve has been outlined.60Transverse Relaxation Time (T2)The transverse relaxation time, T2, represents the gradual reduction of the transversecomponents of the magnetic moment to zero as the nuclei seek equilibrium with the external field,H0. This process occurs through the dephasing of the spinning nuclei relative to each other and isnot accompanied by an energy exchange with the lattice. Thus it is commonly referred to as thespin-spin relaxation. The equation representing the decay of the transverse component ofmagnetization to zero is:f M,y = - ;j!:— [2]The spinning nuclei dephase after the removal of the if field due to variations in themagnetic field experienced by the nuclei. The field experienced by each nucleus is a function of theexternal field and the surroundings of the nuclei. The surroundings of the nuclei vary widely,therefore the field and the precession frequencies of the nuclei also vary.The primary difficulty associated with measuring the exact spin-spin relaxation time is thatthe imposed field, Ho, is usually inhomogeneous and this has a larger effect on the dephasing ofthe spins than does their immediate surroundings. This effect is routinely removed using theCPMG, or the Carr-Purcell-Meiboom-Gill, pulse sequence. This pulse sequence Consists of one90° pulse followed by a sequence of 1800 pulses separated by time t. The initial 900 pulse sets upthe magnetization perpendicular to the applied field and then, due to inhomogeneities in the field,some of the nuclei precess faster, some slower and the dephasing is quite rapid. However, one180° pulse applied at time ton the x axis will reverse the magnetization along the Y axis and haveno effect on the magnetization along the X axis. This pulse moves the magnetic moment of thosenuclei which precess faster to a position behind the average and those which precess slower to aposition ahead of the average. Therefore, after another time interval t, the magnetic moments will61re-align and the dephasing, which was not a result of the inhomogeneous magnetic field, may bemeasured.Spin-spin relaxation times are shorter in liquids confmed in solids; however, the datacollected is a function of the pulse separation time, t, and data sets may only be compared if theeffects of the t variation may be accounted for (Hedberg, 1990, and Kleinberg and Horsfield,1990). In addition, an investigation into the relation between the spin-spin relaxation data and poresizes by Kleinberg and Horsfield (1990) indicates that this relation is more complex than thecorresponding relation between spin-lattice relaxation data and pore sizes. Therefore, T2 data werenot analysed in this study.The NMR Response of Fluids Confined in RocksThe relaxation time of fluids confmed in porous media are much shorter than the timesmeasured for bulk samples of the same fluid. For rocks, this fact was first demonstrated byBrown and Fatt (1956) and has since been confirmed by other researchers. This phenomenon hasalso been documented in such diverse media as animal tissues, biological cells, cheese, and woodlumens. (Bratton et al., 1965, Menon et al., 1987, and Callaghan et al, 1983)Many explanations of the increased relaxation rate/decreased relaxation time have beenpresented. They include:1) the nuclei at the surface may be “bound” to the surface and thus act more like nuclei insolids, which have a shorter relaxation time than those in liquids2) paramagnetic impurities on the pore surface enhance relaxation3) the contrast between the magnetic susceptibility of the fluid and the solid sets up amagnetic gradient within the poresThe relaxation of the nuclei of interest in this study, protons, has been modelled usingdiffusion and exchange mechanisms (Knispel, 1981, Brownstein and Tarr, 1977 and 1979). Thediffusion model is the most widely accepted model and assumes that an increased relaxation rate at62the surface and the natural diffusion of the nuclei account for the increase in the relaxation rate ofnuclei confmed in porous media. Brownstein and Tarr (1979) divided pore sizes into threeregimes: the fast, intermediate, and slow diffusion regimes. In the fast diffusion regime, the poresare small and all of the nuclei diffuse to the surface quickly enough to produce an averagerelaxation time for the pore. In the slow diffusion regime, the nuclei farther from the surface relaxat rates closer to the bulk relaxation rate of the fluid.Slow DiffusionThe theory of nuclear magnetic relaxation of fluids within pores in the slow diffusionregime was developed by Brownstein and Tarr (1979). Pores within the slow diffusion regimemeet the condition that paiD> 10.0, where p is the strength of the surface to relax the protons andis referred to as the the surface sink parameter, a is the radius of the spherical pore, and D is theself diffusion coefficient of the fluid. Theoretically, pores large enough to be within the slowdiffusion regime are associated with an infinite sum of exponential decays with relaxation times Tand amplitudes I which , for a spherical geometry, are given by:212(sinC - Ccos)ln= [3]- sin(2C)]T = a2/D [4]where Cn are the positive roots of:1-cot=paiD [5]63The expression for the observed magnetization M(t,R) at time t for a pore of size R isM(t,R) = M(O) Ie4fTn [6]where M(O) is the initial uniform magnetization and is a function of the quantity of protons presentin the pore, not a function of R. The relaxation times T and the amplitudes I are functions of thesize of the pore, the diffusivity of the fluid, D, the surface sink parameter, p , and the relaxationtime of the bulk fluid, TbTherefore a rock with a distribution of pore sizes, in which the pores are assumed to beapproximately spherical, the T1 or T2 relaxation curve is the integration of a set of infinite sums ofthe form [6] and is given by:y(t)= j c(R) M(t,R) dR [7]where c(R) is the number of spheres of radius R, where e is the signal density and convertsfrom number of spheres to signal intensity.The equations presented above are generally applicable to any type of medium as theyrequire very few assumptions: they may be applied to pores of any size, and modified to accountfor most geometrical shapes. The primary assumptions associated with these equations are (1) thatthe pores are isolated, i.e. very few or no nuclei from one pore diffuse into another pore during therelaxation process, (2) the pores may be approximated by a simple geometric shape such as asphere or cylinder, and (3) the surface sink parameter, p. is constant throughout the medium.64Fast DiffusionThe fast diffusion regime is a special case of the slow diffusion equations presented above.if the pores are small enough that pa/D << 1.0, diffusion will ensure that all of the nuclei reach thepore surface and relax more quickly than they would independently. This results in an averagerelaxation rate throughout such pores so that each pore is associated with one relaxation time andthe expression for the decay curve of a rock simplifies to:y(t)= J s(T)eti1ldT [8]where y(t1) is the data measurement at time t1, and s(T) is the amplitude of the component atrelaxation time T. s(T) is equivalent to M(0) in equation [1] because it is the initial magnetizationassociated with relaxation time T. The assumptions associated with the fast diffusion equationpresented above are slightly more resthctive than those of the slow diffusion equations. Theassumptions are: (1) the pores are small enough so that pR/D << 1, (2) the pores are isolated, and(3) the surface sink parameter is constant throughout the medium. However, when pore sizes areconsidered, the benefit of the fast diffusion regime is that a general surface area to volume ratio ofthe pore , S/Vpore, may be determined. If required, this ratio may then be used to calculate a poresize given an assumption of shape.The one-to-one correlation between pore sizes and relaxation times in the fast diffusionregime allows the direct calculation of pore sizes from spin-spin or spin-lattice decay curves. Theapproach usually taken to this problem is to use the “two fractions in fast exchange” model firstpresented by Zimmerman and Brittin (1957) and further developed by Brownstein and Tarr(1977). This model assumes the pores are within the fast diffusion regime and that the fluid withinthe pore consists of two distinct phases: a bulk phase with a relatively long relaxation time (Tib)65and a surface phase with a substantially shorter relaxation time (Tis). The expression which relatesthe two responses to the measured relaxation time for each pore is:[9]T1 Tib Tis(Munn and Smith, 1987, Gallegos et aL, 1986) where T1, Tib, and Ti are the relaxation times ofthe pore, the bulk fluid, and the surface layer of fluid respectively. fa and fb are the fraction ofbulk and surface water. The modifications of this relation which have been utilized by researchersare all very similar, ranging from Munn and Smith’s (1987) equation:....L_.L 2 10T1 Tib TlsurfrpwhereT1is actually the relaxation time of the surface fluid divided by the thickness of thesurface layer, 8, and rp is a mean pore radius defined as two times the ratio of pore volume tosurface area,2(V/S)pore, to the expression employed by Schmidt et al. (1986):1(1-SSIV) + os/V [11]T1 Tib Tiswhich does not neglect the fraction of surface water as the previous equation, [10], does.The expressions described above are all forms of the general equation for a spherical pore:[12T1 V TibFor example, in equation [11] the equivalent form would be:1_ 6S1 1 1 . _,1 1v—— )+ with P—11 V ls ‘lb ‘lb ‘is ‘lb66This equivalence is direct as p. 6 and T1s are unknowns which may only be estimated from otherdata.EXPERIMENTAL PROCEDURESTwenty-two carbonate samples were selected for use in this study. These samples werechosen to provide an adequate test of the NMR and imaging techniques in this thesis as well as toenable the determination of what properties of the samples affect the NMR response. Thecarbonates were selected to provide a range of pore sizes, with six of the samples having poresgreater than 0.5 mm in length and three of these with very visible (>1mm) pores. Nineteen of thecarbonates, referred to as S3 1- S49, were from a single formation, the Elkton Member of theTurner Valley Formation. This group of carbonates was chosen to ensure that any variation in theNMR results for these samples would be due to variations in the pore network rather thandifferences in the interaction between the protons and their surroundings. The other three samples,referred to as S54, S55 and S77, were also from the Turner Valley Formation. These sampleswere selected because they possessed large, easily visible pores which were not present in thesamples from the Elkton Member.Sample PreparationThe carbonate samples were prepared and the porosity and permeability measurementscompleted by Shell Canada Ltd. These rocks were one inch diameter plugs which had previouslybeen cleaned using methanol extraction followed by drying in a vacuum oven to remove allhydrocarbons. The samples were trimmed to one inch in length and the ends ground smooth andparallel to each other. The one inch length was required in order to ensure that the entire samplewas within the area in the NMR sample chamber in which the magnetic field was constant.67Sample CharacterizationPorosity and PermeabilityThe porosity of the carbonate samples was measured using a dry weight/saturated weightcomparison. The dry weight of the sample was measured after methanol extraction and vacuumdrying to ensure the samples were completely dry. The volume of the plugs was determined bymeasuring the length and diameter of the samples with calipers, and the saturated weight of thesample was measured. This method assumes that the naptha penetrates all of the “effective” orinterconnected porosity in the sample. The equation used to calculate the porosity was:(Ms-Md)“ D(irr2h)where M5 is the mass of the saturated sample, Md is the mass of the dry sample, r is the radius ofthe sample, h is the length of the sample, and D is the density of the naptha.The air permeability of the carbonate samples was measured by determining the pressuredrop across a sample for a constant flow rate of air. Darcy’s Law was used to calculate thepermeability using the equation:— -kpA(h-h1)where Q is the total discharge of fluid (air) per unit timeA is the cross-sectional area of the flow pathL is the length of the flow pathP is the density of the fluid11 is the dynamic fluid viscosity68h2-h1 is the pressure drop across the flow path,k is the permeability of the sample.The porosity and permeability of the samples are summarized in Table 3.1.Minera1ovThe mineralogy of the carbonate samples was determined by a combination of hand sampleand thin section examination. The thin sections were examined twice: once before and once afterthey were stained to aid in the separation of the carbonate constituents.The stain used was made up of a mixture of hydrochloric acid, Alizarin Red S, andpotassium ferricyanide as described by Lewis (1984) and Evamy (1969). Alizarin Red S andpotassium ferricyanide solutions may be applied separately; however, calcite and dolomite cannotbe distinguished without using both. Table 3.2 indicates the stains which result when solutions ofthese two chemicals are used alone and when they are combined.The process of staining required that the thin sections were clean and dry, the chemicalsmixed with deionized water, and the thin sections immersed in the solutions for one to twominutes. The stained thin sections were dried using compressed air to ensure that the stain did notbecome streaky, that a residue of the stain was not left on the section, and that the section driedrapidly. The use of potassium ferricyanide requires several precautions as it is photo sensitive andwill react with tap water. Therefore only deionized water was used and, if the solution was notgoing to be used within 48 hours, it was stored in a dark bottle away from the light.69TABLE 3.1 Porosity and PermeabilityS3 1S32S33S34S35S36S37S38S39S40S41S42S43S44S45S46S47S48S49S54S55S770.0380.0540.1070.0760.0940.0680.0410.0280.0550.0530.0190.0630.06 10.0600.0590.0490.0660.0490.0420.0570.0480.0820.0500.0700.5500.4100.1300.0600.0400.0500.0600.0500.0400.0500.0500.0500.0500.0500.0500.0500.0500.0700.0605.43SAMPLE POROSiTY PERMEABILITYmD70TABLE3.2StainsforDifferentiatingBetweenCalciteandDolornite*STAININGCALCITEDOLOMITESILICAREAGENTSCompositionsaregiveninFe++freeFepoorFerichFe++freeFe++<1Fe>1weightspercent.CnticalMgMgQUI{’1solutionstrengthsareCHALCEDONYunderlined.CALCITEDOLOMITEFERROANANKERITEOPALFERROANCALCITEDOLOMITEsensustictOsensustricto0.2%hydrochloricacidRedRedRedNotNotNotNotStained0.2%AlizarinRedSStainedStainedStainedQ,2%hvdrochloricacidNotLightDarkNotLightDarkNotStained0.2%poiassiumferricyanideStainedBlueBlueStainedBlueBlue0.2%hydrochloricacidRedMauveVioletNotLightDarkNotStained0.2%AlizarinRedS0.2%potassiumferricyanideStainedBlueBlue *afterEvamy,1969Samples S31 to S49These nineteen samples are from the Mississippian Elkton Member of the Turner ValleyFormation. This formation has previously been described as consisting of three sub-members: alower member of porous, coarsely crystalline dolomite which grades into less porous limestone insome areas, a middle member of dolomitic shale or argillaceous and silty dolomite, and an uppermember of finely to coarsely crystalline dolomite containing varying amounts of chert (Penner,D.G., 1957). This indicates that there is a wide variation in the mineralogy and physical propertiesof rocks from this formation; however, from the literature it would be expected that the rockswould be at least partially dolomitized, the fossils present will be dominated by crinoids and thatother fossils such as bryozoans, brachiopods, and relict texture suggesting pelletoids may bepresent (Penner, D.G,. 1957, Young, F.G., 1967, and Thomas and Glaister, 1960).The hand sample and thin section examinations revealed the samples were primarily matrixsupported limestones which were partially dolomitized, the dolomite being present as small rhombs(<0.2 mm) within the matrix and in partially replaced fossils (Figure 3.2). A small component ofshale (<10%) was also present and the samples were gray, with little to no visible porosity. Whatlittle porosity was visible was present in samples S31, S32 and S33. The fossils observedconsisted primarily of crinoid ossicles and stems (Figure 3.3) with a small number of bryozoansand brachiopods (Figure 3.4). Throughout the nineteen samples, only one area with a relict textureof pelletoids was viewed. A general description of these samples is: grey, low porosity, partiallydolomitized, matrix supported crinoid limestones.Samples S54, S55, and S77These three limestones were from the Mississippian Turner Valley Formation. They werecream coloured and contained a great deal of visible porosity. Pores up to 1 mm in diameter weremeasured in samples S54 and S55 and 2 to 3 mm in diameter in sample S77. The samples were72Figure 3.2: Partial dolomitization in Ellcton Member samples, S31 to S49. Dolomite is blue andusually rhombic, calcite is red or pink, and extremely small crystals consisting of carbonate andsilicate minerals appear as dark areas.73• ... ... --.. .1’4 . •••: •.-.7A.•0•-I.--S0.0 0.25mm•••••-.;,• -.•..•• • •.1..;*Figure 3.3: Examples of crinoid stems and ossicles in Elkton Member samples, S3 1 to S49.Dolomite is blue and usually rhombic, calcite is red or pink, and extremely small crystalsconsisting of carbonate and silicate minerals appear as dark areas. Lower photo contains a crosssection of a crinoid stem; examples of crinoid particles in the upper photo are indicated by C ‘s.74Figure 3.4: Examples of bryozoan and brachiopods which occur in Elkton Member samples, S3 1to S49. Dolomite is blue and usually rhombic, calcite is red or pink, and extremely small crystalsconsisting of carbonate and silicate minerals appear as dark areas.75primarily coarsely crystalline limestone and contained less than one percent silica. S54 and S55contained little or no dolomite, while S77 contained approximately 10% dolomite.DATA COLLECTIONSample PreparationFully Saturated SamplesThe samples were dried by heating them in a vacuum oven for several hours to remove thenaptha after the porosity was measured. The samples were then immersed in 100,000 ppm NaClfor several hours at room temperature and pressure, a vacuum was applied for several hours, andfinally, the samples and fluid were transferred to a pressure vessel and 600 psi of Nitrogen wasapplied for approximately fifteen minutes. Following saturation with the brine, the samples wereimmersed in 100,000 ppm NaCl within a sealed dessicator for storage. The dessicator was kept ina heat bath which was maintained at 40°C for at least one-half hour prior to any NMRmeasurements. The samples were heated so they would be in equilibrium with the NMR samplechamber which is maintained at 40°C. Spin-lattice relaxation times vary with temperature, soitwas important to ensure the temperature of the sample was constant throughout the experiment.The sample to be measured was removed from the dessicator using tongs, surface dried toremove excess water, and then wrapped in 2 layers of cling film. Teflon tape was wrapped aroundthe outside of the cling film in two alternating layers. The layering of cling film and teflon tapewas used to ensure that a minimum of fluid loss from the sample during the measurements. Ateflon filler surrounded the samples within the NIVIR sample holder. The teflon filler, which doesnot produce an NMR signal, was required to keep the samples centered within the holder whichwas approximately 3.8 centimeters in diameter, 1.27 centimeters larger than the samples. The76sample holder was a glass, square-ended test tube which was placed within the sample chamber ofthe NMR machine to collect both T1 and T2 decay curves.Partially Saturated SamplesSix of the samples employed in the experiment on the fully saturated samples were cleanedto remove the brine and re-saturated with deionized water. The brine was removed by soaking thesamples in deionized water for a period of several days, with the deionized water being replacedwith fresh water daily. The saturation with de-ionized water was accomplished in the same waythey were previously saturated with brine.The samples were saturated with deionized water because the variation in saturation wasaccomplished using evaporative drying, which would cause an increase in the salinity of the porefluid if a brine were used. Evaporative drying was accomplished by placing samples in adessicator with anhydrous CaSO4after the initial NMR measurements of the samples in a fullysaturated state. Prior to the fully saturated measurement, the samples were stored in a sealedvacuum dessicator within a heat bath which was turned on for at least thirty minutes prior to themeasurement of the NMR decay curves of a sample.The partially saturated samples were weighed before and after each NMR measurement todetermine their saturation level. The preparation for the NMR measurement was the same as forthe fully saturated samples: if necessary, the sample was surface dried, and then two layers ofcling wrap and two of teflon tape were wrapped around the sample. The cling wrap and teflon tapereduced fluid loss from the samples both during the NMR measurement and as they were heated to40°C. The temperature of the sample was increased to 40°C by placing it in a sealed NMR sampleholder and placing the holder in a 40°C heat bath for at least twenty minutes prior to the NMRmeasurements.77NMR MeasurementsThe NMR measurements were performed with a Bruker Minispec PC-20 NMRSpectrometer with a magnetic field strength of 0.47 T. The pulse sequences and parameters usedto collect each data set were determined by two different EDMs (Experiment Definition Modules)and several calibration tests were performed regularly to ensure the NMR machine was operatingproperly. The homogeneity of the magnetic field, the diode offset and the phase of the phasesensitive detector were checked at least once a week. To ensure the highest possible signal to noiseratio, the built-in bandwidth filter was set to the lowest possible value of 150 kHz.Spin-Lattice Decay CurvesThe spin-lattice relaxation curves were collected using a Bruker EDM 510 which wasprogrammed with the “inversion recovery” sequence described previously. A maximum of 160measurement times were possible, with the time of each measurement calculated with the equation:t=CT1*CFCF varies from 0.04 to 5.12 and is determined by multiplying the previous value by:(n-i)— Vo.oi.where n is the number of points on the curve and the initial value of CF is 0.04. Although theexact time calculated in this manner is presented as the measurement time, the actual measurementtimes of the Minispec PC-20 occur at whole number millisecond times, thus if the Minispec claimsto have measured the magnetization at 2.2 or 2.8 ms, the measurement was actually performed at 2ms. CT1 is a constant which is an estimate of the average T1 of the sample, thus the time between78two “inversion recovery” pulse sequences is set to 5 times this value to ensure that the nuclei havetime to relax completely between measurements. Between 3 and 45 measurements were taken ateach time and stacked to produce the final result. The measurements were completed inapproximately one hour for all samples since for high porosity samples, the data were clearer andthe same data quality was obtained from a smaller number of measurements.Spin-Spin Decay CurvesExperiment Defmition Modules 610 and 612 were used to collect the T2 decay curves ofthese samples. Both EDM’s employ the CPMG pulse sequence described in the theory, however,610 measures the magnetic field strength after every tenth echo while 612 measures after everysecond echo. Thus, 610 is suited to samples with long T2s and 612 to samples with short T2s.Unfortunately, only 50 measurements of the curve may be taken with the 612 EDM while EDM610 will collect up to 160 measurements.The nuclei must be allowed to reach equilibrium with the external field between each spin-echo experiment, thus the constant mentioned above, CT1, is used to ensure that 5 times theapproximate spin-lattice relaxation time separates each spin-echo experiment. The time over whichthe curve is measured is set with the CT2 constant, which should be approximately equal to theactual T2 of the sample. The maximum time at which data is collected is greater than 2*CT2,ensuring that 90% of the decay curve with a time constant of CT2 would be collected. .CT2 wasadjusted to collect the curve at its closest approach to zero; however, it was maintained at a highenough value so that noise effects did not cause the signal to decrease below zero volts.DATA ANALYSISThe spin-lattice relaxation data collected for this study were analysed to determine theporosity, pore sizes, and T1 distributions of the samples. Pore size distributions were calculated79using two methods: linear inversion of the relaxation data using the slow diffusion model, anddirect calculation from the T1 distributions using a relation, equation [12], developed using the“two fractions in fast exchange” model.Data InversionSpin-Lattice Relaxation Time (T1)DisthbutionsThe linear inversion methods presented by Whittall and MacKay (1989) were used todetermine both continuous and discrete T1 distributions. These inversions assumed, as is the casein either fast or slow diffusion, that the decay curve produced by a fluid saturated rock will bemultiexponential in form. Thus, the equation used to determine the appropriate distribution of Ti’swas equation [8]:y(t1)= J s(T)etiul’dT [13]awhere y(t) is the datum at measurement time t1, and s(T) is the amplitude of the component atrelaxation time T.The linear inversion method assumes a number of known relaxation times Tj and solves forthe associated amplitudes s(T) or sj. In order to execute the inversion with a computer, theequation above was discretized assuming the model was a sum of M delta functions with areas sj atknown relaxation times T and the resulting expressions are:s(T) = s(T-T) [14]80My1=sjgjj i=1,...,Nj= 1where the kernel, gjj, is given by:fTj+lgjj=etiIJdTJTThe non-negative least squares (NNLS) algorithm selected to invert the data would providethe closest possible fit to the data by calculating the solution which minimized the term:E I sjgjj -i=1 j=1However, experimental data always contains a certain amount of noise, therefore the preferredsolution should not fit the data exactly. Solutions to the inversion should misfit the data by anamount corresponding to the error in the data, which is represented by the standard deviation, a.2This is accomplished usmg the X statisticN= (y - y)Ia [15]to quantify the misfit between the measured data, y, and the data predicted by the constructedmodel, y. The expected value of is equal to the number of data points, N, corresponding to amisfit of approximately one standard deviation for each data point.The NNLS algorithm employed the criteria as a measure of the appropriateness of themodel, s(T), produced. Two approaches were implemented: one which assumed that theamplitudes were delta functions and produced a discrete spectrum as described above, and one81which incorporated additional constraints in order to produce a piece-wise continuous distribution.The piece-wise continuous distribution was obtained by calculating the solution which minimizedthe term:sjgjjyj2+p.s?i=1 j=1 j=1in which u is a trade-off parameter which, if 0, allows the NNLS solution to be calculated, and asit becomes larger, places more emphasis on a solution which minimizes the “energy” of thesystem. Minimizing the second term, referred to as the energy of the system, corresponds tominimizing the structure in the solution.Pore Sizes - Slow Diffusion RegimeThe NNLS inversion techniques described above were also employed to determine poresizes based on the theory of Brownstein and Tarr (1979). The pores were assumed to be sphericalin shape and the inversion equation, [7], was:y(t) = j c(R) M(t,R) dRThe application of this equation was accomplished by discretizing the problem, providing aset range and number of pore radii as a basis for the inversion, and removing the constants fromthe equation and solving a simplified equation:y(t) = f c(R) R3 M(t,R) dRJo82These radii and equations [3], [4], and [5] were used to calculate the distribution ofamplitudes and relaxation times associated with the selected pore sizes. These amplitude/relaxationtime distributions were then utilized to produce a best fit model to the NMR decay curve bydetermining the number of pores of each pore size present in the sample.The development and testing of the computer algorithm which implemented this inversionis described in further detail by Whittall (1991).Porosity DeterminationThe porosity of the samples were determined based on a basic principle of NMR: themagnitude of the initial magnetization, M(O), is proportional to the number of protons present inthe sample. Therefore, if a sample is fully saturated, the “effective” porosity can be determined bycomparing the initial response of the sample to that of a fluid sample of the same volume.Similarly, the saturation of a partially saturated sample can be determined by comparing itsresponse with that of the fully saturated state.The magnitude of M(O) can be determined from the T1 distributions of the samples bysubstituting t =0 into equation [13], producing the expression:y(O)=f s(T)dT [16]JoThis indicates that the magnitude of M(0)=y(O), and thus the relative quantity of fluid, may becalculated by integrating over the T1 distribution determined for a fully or partially saturated rocksample.83Pore Sizes - Fast Diffusion RegimeThe pores sizes of the samples were calculated using the “two fractions in fast exchange”model by using equation [12]_LT1 V T,Throughout the remainder of this thesis, the calculation of pore sizes from the T1 distributionutilizing this equation will be referred to as the fast diffusion method.RESULTSFully Saturated SamplesI1DistributionsContinuous and discrete T1 distributions were produced for each fully saturated sampleusing the non-negative least squares inversion method described earlier. Graphs of thesedistributions are included in Appendix 6.The discrete relaxation spectra for samples S34 to S49 are very similar. As illustrated bythe example shown in Figure 3.5a, most of the fluid in the pores had relaxation times just under 10ms and the maximum relaxation time ranged from 60 to 200 ms. The inversion of data whichproduced a non-zero amplitude at relaxation times greater than 90 ms was difficult. Regardless ofthe range of relaxation times chosen to invert the data, a peak was always present at the maximumtime chosen.84Samples S3l to S33 were slightly different than the others from the Elkton Member, asshown by the S3 1 spectra in Figure 3.6a. The maximum amount of fluid in the pores had arelaxation time between 40 and 60 ms and, for S31 and S32, a substantial amount of the fluid(>20%) possessed a relaxation time of 300 and 200 ms respectively.The continuous relaxation spectra for the samples from the Elkton Member are exemplefiedby the S39 spectra in Figure 3.5b. These spectra all reached a maximum at 20 ms or less;however S31,S32 and S33 were represented by slightly larger T1 values and attained theirmaximum amplitude between 25 and 80 ms. The increase in relaxation times for S31, S32 andS33 relative to the other Elkton Member samples is illustrated by the continuous relaxation spectraof S31 shown in Figure 3.6b.The discrete and continuous spectra of S54,S55 and S77, of which the S77 spectra inFigure 3.7a and 3.7b are typical, are represented by relaxation times close to a full second. Thediscrete spectra were characterized by maximum amplitudes at the maximum relaxation timesrepresented: 1750, 1000, and 1250 ms for S54, S55 and S77 respectively. The continuousspectra for these samples were similarly biased towards the longer relaxation times, with themaximum in the T1 curve occurring at 800, 750 and 1000 ms for S54, S55 and S77 respectively.PorosityThe total integral of the T1 spectra, M(0), of each sample should be related to the porosity of thatsample, as shown by equation [15]. This correlation is shown in Figure 3.8 for both theContinuous and discrete relaxation time spectra. The best fit straight line for the graph of porosityvs M(0) has a regression parameter of 0.728 for the continuous T1 spectra and 0.883 for thediscrete spectra.85Discrete Relaxation SpectrumSpin-Lattice Relaxation Times (ms)0.60.50.40.30.20.10.0=1 10 100 1000It)Spin-Lattice Relaxation Times (ms)(a)Continuous Relaxation Spectrum0.080.060.040.020.00 ...I1 10 100 1000(b)Figure 3.5: Sample S39 Ti distributions.(a) discrete and (b) continuous.86Discrete Relaxation Spectrum0.4j0.1•0.01 10 100 1000Spin-Lattice Relaxation Times (ms)(a)Continuous Relaxation Spectrum0.081 10 100 1000Spin-Lattice Relaxation Times (ms)(b)Figure 3.6: Sample S31 Ti spectra.(a) discrete and (b) continuous.87Discrete Relaxation Spectrum(a)Continuous Relaxation Spectrum0.0810 100 1000 10000(b)Figure 3.7: Sample S77 Ti spectra. (a) discrete and (b) continuous, interpolate a coefficient of2.25x i0- m2/s at 40°C, the temperature of the sample chamber in the NMR spectrometer used tocollect the data for this study.a)a)a)0.80.60.40.20.01 10 100 1000 10000Spin-Lattice Relaxation Times (ms)0.06’0.04’0.020.00 ‘I ...ISpin-Lattice Relaxation Times (ms)88Correlation of Porosity with M(O)Continuous Ti Spectra12000 y= .-163.36+9.4015e+4x RA2=07288000M(O) 60002000 EEEEEZEZEZ40000.00 0.02 0.04 0.06 0.08 0.10 0.12PorosityCorrelation of Porosity with M(O)Discrete Ti Spectra12000- y = 834.74 + 1.0175e+5x R”2 0 883100008000M(O) 6000400020000- • • I • I • I •0.00 0.02 0.04 0.06 0.08 0.10 0.12PorosityFigure 3.8: Coffelation of the integral over the the discrete and continuous Ti distributions(M(0)) with measured porosity.89Pore SizesThe pore sizes calculated using both the slow and fast diffusion methods require theprior knowledge or estimation of three parameters: D, the self-diffusion coefficient of the bulkliquid; p, the surface sink parameter; and Tib, the relaxation time of the bulk fluid. The surfacesink parameter, p, may only be determined through comparison of the calculated pore size spectrawith a previously determined pore size distribution; however, both D and Tib may by measured.The relaxation spectrum for the 100000 ppm NaC1 brine used to saturate all of the carbonatesamples is shown in Figure 3.9. Tib was determined to be 2.0 s for water which had not been degassed to remove the paramagnetic oxygen, a value which compares favourably to othermeasurements (Schmidt et al., 1986, Endom et al., 1967). The self-diffusion coefficient may onlybe measured with a magnetic gradient, which was not available, therefore the values determinedby Endom et al. (1967) for water molecules in 100,000 ppm NaC1 at 0, 25, 50 and 80°C were usedto interpolate a coefficient of 2.25xiO m2/s at 40°C.Fast DiffusionThe fast diffusion method was applied to calculate those pore radii which corresponded tothe relaxation times determined with the NNLS inversion (Appendix 6). The surface sinkparameter, p. was calculated by determining the value required to equate the average T1 of thediscrete and continuous S55 and S77 spectra with the average pore size determined with theconfocal laser scanning microscope. The values calculated were lx10 m/s for three of thecalculations, and 2x105mis in one case, therefore the most common value, lxl&5mis, wasselected. When the T1 distributions were altered to pore size distributions, the discrete spectra wereneglected because the porosity in most rocks is better represented by a continuous distribution.90Bulk Relaxation Time of 100000 ppm NaC11.0•0.841.)0.2 -0.0- •__.1 10Spin-Lattice Relaxation Time - Ti (s)Figure 3.9: 100,000 ppm NaC1 brine relaxation time.91The two pore size ranges determined for samples S3 1 to S49 using the fast diffusion poresize calculations and the relation between the original T1 spectra and the final pore size spectra areillustrated by the T1 and pore size spectra of S39 and the pore size spectra of S31 shown in Figure3.10. The most striking result for samples S31 to S49 was the direct alteration of relaxation timesto pore sizes with no change in the overall shape of the distribution. The maximum amplitude ofthe disthbutions covered two size ranges: between 0.3 and 0.7 p.m for S33-S49 , and 1.0 to3.0 p.m for S31 and S32.The T1 distributions of S54, S55 and S77 contained much larger relaxation times, thereforethe resulting pore size distribution was “stretched” as the relaxation times in the original distributionbegan to approach the relaxation time of the bulk fluid (Figure 3.11). The maximum amplitude ofthese spectra occurred between 40 and 60 p.m. sizes ten to twenty times larger than the radii atwhich the maximum amplitudes of S31 and S32 occurred.Slow DiffusionThe ability of the slow-diffusion inversion method to reproduce the correct pore sizedistributions was tested by producing a theoretical spin-lattice decay curve based on a theoreticalpore size distribution similar to that of sample S77. The times at which the theoretical“measurements” were calculated were varied by matching the actual measurement times of thesample S77 data ( 38 to 4863 ms) and by expanding this to include a much larger range ofmeasurement times, from 1 to 22000 ms. The inversion of the theoretical data matched thetheoretical distribution more accurately when the time range of the data was expanded, as illustratedin Figure 3.12.The theoretical data was also altered by adding 1% and 5 % noise and the subsequentinversion produced pore size spectra which matched the actual distribution reasonably well. Figure3.13 illustrates the loss of detail and accuracy with increasing noise through comparison plots of92Ti and Pore Size Distributions0.08Pore Sizes Ti0.06.0.040.020.00• •..—.. ..‘ —I I ••• ...i08 io 10.6 i05 i0 i0 10.2 10i io°Ti (ms) and Pore Radius (m)(a)Pore Size SpectraFast Diffusion Method0.080.060.00108 106Pore Radius (m)(b)Figure 3.10: (a) Sample S39 Ti and fast diffusion pore size distributions. (b) S31 fastdiffusion pore size distribution.93Ti and Fast Diffusion Pore Size Spectra0.080 Pore SizesI.Ti0.04’I0.02o.oo —.i0 i0 io io’ io i02 10’ 100 iOiTi (s) and Pore Radius (m)(a)Pore Size SpectraFast Diffusion Method0.060.05C.—0.04tl) 0.03C0.01W0.00•106 i-Pore Radius (m)(b)Figure 3.11: (a) Sample S77 Ti and fast diffusion pore size distributions. (b) S55 fastdiffusion pore size distribution.94Slow Diffusion and Confocal Pore Size Spectra0.4 Theoretical Pore Size SpectrumSlow Diffusion Pore Size SpectrumTheoretical data range: 38-4863 ms0.30.—I. I0.2 $ ik_I’ IIo IjIi I0.1_________iLl .io-3Pore Radius (m)Slow Diffusion and Confocal Pore Size Spectra0.25Theoretical Pore Size Spectra0.20 Slow Diffusion Pore Size SpectraTheoretical data range: 1-22000 ms..II0.15- :II:: :0.050.00 - — —f •—i—’.! 1io5 io ioPore Radius (m)Figure 3.12: Slow diffusion method results of inverting theoretical data produced using atheoretical pore size spectrum similar to that of S77. Variation in the pore sizes is a result of usingtwo time ranges to calculate the theoretical decay curve. The upper graph was produced with dataover the same time range as the data for S77 was collected, the lower with an extended range.95Slow Diffusion and Confocal Pore Size Spectra0.40.3_____Slow Diffusion Pore Size SpectrumTheoretical data range: 38-4863 msSlow Diffusion Pore Size Spectrum0.2 Theoretical data range: 38-4863msNoise: 1%II I________ __io io3Pore Radius (m)Slow Diffusion and Confocal Pore Size Spectra0.6Slow Diffusion Pore Size SpectrumTheoretical data range: 3 8-4863 ms0.4 II Slow Diffusion Pore Size SpectrumTheoretical data range: 38-4863msII Noise: 5%• 0.3S 0.20.1 Ri i______ ____oo —- —io4Pore Radius (m)Figure 3.13: Slow diffusion method results of inverting theoretical data based on a theoreticalpore size spectrum similar to that of S77. The pore sizes resulting when 1 and 5% noise areadded to the theoretical data are shown in the upper and lower graphs respectively.96the pore size spectra calculated from the perfect and noisy theoretical data. The pore size spectradetermined from noisy data are characterized by the appearance or enhancement of the amplitudeassociated with the largest pore size and the representation of the pore size spectra by two or threepore sizes rather than a continuum of pore sizes. If the component at the largest pore size wasincluded, the average pore size determined from the data with 1% noise was 153 pm, and from thedata with 5% noise was 216 pm. This illustrates the sensitivity of the problem to additional noise,as 153 .tm is within 20% of the actual average of 122 p.m, while 216 p.m is almost double thisexpected value. The results when the final pore size is discounted are slightly more promising,with average pore sizes of 106.1 and 82.3 p.m for the data with 1% and 5% noise respectively.The variation in the X2 parameter with p was tested with the theoretical data. Whittall2(1991) found that the minimum mX tends to underestimate the actual parameter, particularly forpore sizes in the slow diffusion regime. For the theoretical data produced for this study, theinversionof theoretical data without added noise produced a minimum X2 when p matched the value the datawas produced with. When even 1% noise was added to this data, this minima was either alteredfrom an obvious to a broad minimum, or, alternatively, the value decreased steadily as p wasdecreased.When applied to both the real and theoretical data, the slow diffusion method consistentlyproduced spectra which more closely resembled discrete than continuous pore size distributions.Typical results for the theoretical data are shown in Figures 3.13, while examples of the spectracalculated from real data are shown in Figures 3.14 and 3.15. The surface sink parameter used foreach real data set was varied and the situations described for the theoretical data, in which a broadminimum in or a constant decrease in with p occured, were also observed for the real data.In all cases, either the p value just above the minimum, a p value slightly smaller than that at whichincreased substantially or, if corroborating data was available, the p value which best fit theknown pore size distribution was selected. The surface sink parameter for samples S3 1 to S49ranged from 1.5 to 3.75x 10-5 mIs, with an average of 3.0 x105 ni/s. The minimum occured at97p values of 2x105 rn/s for S54 and S55, and 3.5x10 rn/s for S77. However, surface sinkparameters of lx iO rn/s and 1x10 rn/s fit the confocal pore size distribution more accurately forS55 and S77 respectively and therefore these values were used for these two samples.The pore size spectra produced for all of the data sets possess a peak at the largest pore sizeincluded in the inversion range. As shown in S31 and S43 pore size spectra in Figures 3.14 andthe S55 and S77 spectra in Figure 3.15, this peak represents from 15 to 25 percent of the integralover the complete spectrum for samples S3 1 to S49 and between 50 and 60 percent of the completespectra for samples S54, S55 and S77.The pore size associated with the maximum amplitude of the S33 to S49 distributions ranged from0.8 to 1.5 microns with several samples having a second peak at pore sizes between 1.5 and 6microns, as illustrated by sample S43’s pore size spectra in Figure 3.14a. The spectra producedfor S3 1 and S32 were represented by larger pores than the other samples from the Elkton Member.This increase in pore sizes relative to the S33 to S49 results is typified by the pore size spectra ofS31 in Figure 3.14; the maximum amplitude and secondary peaks occur at 2.5 and 10 p.m for S31and at 4 and 15 p.m for S32.Two peaks other than the one present at the largest pore size were produced in thedistribution of each of the other three samples, occurring at 3 and 9 p.m for S54, at 1,5.2, 6.5, 37,and 43 J.tm for S55, and at 10 ,50, and 59 p.m for S77. These results are illustrated by the poresize spectra of S55 and S77 included in Figure 3.15.The pore size spectra described above are non-unique solutions to a complex problem, and,since the pore size and amplitude of the distribution are affected by the parameters required for theinversion, the variation in the spectra with D, p and Tib was determined. Several sets of data fromall of the samples were tested and samples S3 1 to S49 were found to have similar responses, asdid samples S54, S55 and S77. An example of the variation in the pore size spectra for the twotypes of responses to Tib is shown in Figure 3.16. The spin-lattice relaxation time has only a verysmall effect on the amplitude of the spectra for S3 1 to S49, but it causes an increase in pore sizes98Pore Size SpectrumSlow Diffusion Method0.300.2I0.0 •, . A.,i06 io5Pore Radius (m)(a)Pore Size SpectrumSlow Diffusion Method0.—0.2I0.1•0.0• IIi0 106 io5 1oPore Radius (m)(b)Figure 3.14: Slow diffusion pore size distributions. (a) S43 (b) S3199Pore Size SpectrumSlow Diffusion Method0.60.5C0.4j :::::z-iA..106 ioPore Radius (m)(a)Pore Size SpectrumSlow Diffusion Method0.8C.—0a) 0.40•1Pore Radius (m)(b)Figure 3.15: Slow diffusion pore size distributions. (a) S55 (b) S771000C)0Figure 3.16: Examples of the variation in the pore sizes determined using the slow diffusionmethod when the value of the bulk fluid relaxation time, Tib, is varied.as it is decreased forsamples S54, S55 and S77.Pore Size Variation with the Bulk Relaxation TimeSample S3920I.1•4.Os•2.Os-1.0s0•0.48 1.69I1.98Pore Size0.80.62.31 45.63Pore Radius (m)Variation with the Bulk Relaxation TimeSample S554.Os2.Os•0.4020.0-Pore Radius (m)101The effects of varying D are illustrated in Figure 3.17 by the variation in pore sizes forsamples S39 and S55. For samples S31 to S49, the variation with D was not direct since the poresizes did not change according to a pattern, and the overall average pore size remained relativelystable. However, as D was increased for samples S54, S55 and S77, the pore sizes also increased.The surface sink parameter has a defmite effect on the pore size spectra determined for all of thesamples; as shown in Figure 3.18, a larger surface sink parameter produces a larger pore size. Theamplitude of the peaks vary as well, with the amplitude at larger pore sizes increasing with thesurface sink parameter.Partially Saturated SamplesFluid DistributionsThe Ti spectra corresponding to varying saturations of deionized water for six samples(S31,S41,S42,S44, S54 and S55) were determined and are included in Appendix 7. As shown inthe S54 spectra in Figure 3.19, the T1 spectra of the six samples generally decreased in amplitudeand shifted to the left towards smaller T1s as the saturation of the samples was decreased. Theshift towards higher T1swhich was observed in varying degrees at the lowest saturations ofsamples S41, S42, S44, and S55, is illustrated by the T1 spectra of sample S44 in Figure 3.20.Fluid SaturationsThe integral over the T1 distributions (M(0)) was correlated with the saturation of thesamples and, as shown in Figure 3.21 for samples S31 and S41, the relationship could beapproximated by a straight line. Graphs of M(0) vs saturation for all six samples are in Appendix7 and the regression parameter indicating how well the data fit a straight line ranged from 0.92 1 to0.993.1 02C0a.)EC>Pore Radius (m)Pore Size Variation with the Self-Diffusion CoefficientPore Radius (m)00 C’ oOCr C— — C4 \OFigure 3.17: Examples of the variation in the pore sizes determined using the slow diffusionmethod when the value of the self diffusion coefficent of the bulk fluid, D, is varied.1 .21.0 -I0.8-I0.6-IPore Size Variation with the Self-Diffusion CoefficientSample S39J 4.5e-9m/s2.25e-9 m2/sR 1.125e-9m/sI___0.4 -,0.2-’0.48 1.69 1.98 2.31 45.63Sample S55C0a)—C1.51.0:: I—F4.5e-9m21s2.25e-9 mI 1.125e-9m/s1 03Pore Size Variation with the Surface Sink Parameter0j00E01.0Sample S390.0 -,00 .:c cn c’iRadius (m)7.5e-5 rn/s3.75e-5 rn/s•00 \C to o oPorePore Size Variation with the Surface Sink ParameterSample S5524e-4 rn/s1- 2e-4 rn/sle-4 rn/s0- ..-----.—--. • —OcCC00O00Pore Radius (m)Figure 3.18: Examples of the variation in the pore sizes determined using the slow diffusionmethod when the value of the surface sink parameter, p, is varied.1 04SampleS54NormalizedTiSpectraatPartialSaturationsSaturations0.995a0.920—0-——0.886—1.———0.7940.730—0--—0.6320.5570.4640.362-1-0.2370 (7’Spin-LatticeRelaxationTime(ms)Figure3.19:Exampleof thedecreaseinamplitudeandrelaxationtimewithsaturationforpartiallysaturatedsamples0.08 0.060.040.020.0010100100010000SampleS44NormalizedTiSpectraatPartialSaturations0.060.05Saturation0.040.9140.8120.03---0.711—0——0.6140.020.509—0——-0.3240.01 0.001000Spin-LatticeRelaxationTime(ms)Figure3.20:ExampleoftheincreaseinlongT1valuesatlowsaturationsrelativetoT1values athighersaturations.1101000)Sample S31Correlation of Saturation with M(O)22002000y—41012+21238x R”2_09211800M(O)160014001200 •0.6 0.7 0.8 0.9 1.0SaturationSample S41Correlation of Saturation with M(O)1300- y = 14.402 + 1203.lx R’2 = 0.99312001100M(O)600’’9008007000.5 0.6 0.7 0.8 0.9 1.0SaturationFigure 3.21: Examples of the correlation between fluid saturation and the total integral over the Tidistribution (M(0)).1 07DISCUSSION AND CONCLUSIONSThe T1 distributions determined for porous samples are expected to reflect the pore sizes ofthe samples, with larger T1 values indicating the presence of larger pores. This expectation wasmet for both the discrete and continuous T1 spectra. Considering only those samples from theEllcton Member, the porosity of these samples was primarily microscopic porosity, with S31, S32,and S33 possessing minor amounts of visible porosity. Correspondingly, the continuous anddiscrete T1 spectra of samples S43 to S49 reached their maximum amplitudes at relaxation times of12 ms or less, while samples S31 to S33 reached their maximum amplitudes at relaxation timesbetween 25 and 80 ms. The correlation between pore sizes and relaxation times was furthersupported by the continuous and discrete T1 spectra of samples S54, S55 and S77 possessingmaximum amplitudes at 750 to 1750 ms. The discrepancy between the magnitude of the S31 toS49 relaxation spectra and those of samples S54, S55 and S77 was expected because the majorityof the porosity of the final three samples was visible. The average pore sizes determined forsamples S55 and S77 using the confocal imaging technique discussed earlier were 25 and 122 jim,which compare qualitatively with the average relaxation times the discrete and continuous T1spectra occurring at 580 (discrete)! 470 ms (continuous) for sample S55 and 909/1392 ms forsample S77The correlation between porosity and the integral over the T1 distributions (M(0)) was notvery accurate, with regression parameters for the best fit straight lines of 0.728 and 0.883 for thecontinuous and discrete spectra respectively. There are three possible causes for this disparity inthe theoretical and experimental correlations between porosity and M(0). Firstly, the samples areassumed to be fully saturated when in reality they are not, and furthermore, each sample will not besaturated to the same extent, therefore a certain amount of scatter is expected. The variation insaturations and the expectation that the samples are not fully saturated was formed by measuringthe dry and saturated weights of several of the samples, calculating the volume of fluid in the108sample, and comparing this to the expected volume of the porosity. In most cases, the porosity inthe samples actual saturation ranged between 80 and 90%. The second possible source of error inthe correlation between porosity and M(0) lies in the accuracy of the porosity measurements: withthe method used, error in the porosity may be as large as 0.5% porosity. The final contribution toerror may be attributed to the difficulty in determining the actual value of M(O). Many of the decaycurves are very steep at short measurement times, thus small changes in the fit of the inversion datato this portion of the curve can result in substantially different values of M(O).The pore sizes calculated using both the slow and fast diffusion methods agreedqualitatively with the presence and/or absence of visible porosity in the samples, with pore sizesfrom smallest to largest being calculated for the samples in this order: S34 to S49, S33, S31, S32,S54 and S55, and S77. Quantitatively, the comparison between both fast and slow diffusion poresizes and the confocal pore size spectra shows reasonable agreement between the three. Thiscomparison is shown in Figure 3.22 for samples S55 and S77. The inaccuracy of the fastdiffusion pore size spectra lies primarily in the stretching of the distribution at the highest poresizes due to the fact that the equation used to calculate pore sizes from T1 values becomesundefmed as T1 approaches Tib. This occurs due to the term:(1 1 ‘i-i‘T1 Tib’in the equation used to calculate pore sizes from T1 values with the fast diffusion method.The primary difficulty involved in interpreting the inversion results for the samples lies inthe determination of the surface sink parameter, p. The values of p used in the fast diffusion poresize calculations were an order of magnitude less than those which produced the best fit using theslow diffusion method, and, as indicated by the inversion of theoretical data with noise added to it,determining the correct p value is very difficult unless an independently determined pore sizedistribution is available. The problem just mentioned may, if the data are accurate, be turned intoan alvantage since a minimum in will occur slightly below the actual value for the sample and1 09thus the surface sink parameter can be determined without the availability of an independent poresize disthbution. Thus the slow diffusion method can be more useful than the fast diffusionmethod.One of the objectives of this study was to determine whether the fast diffusion method wasadequate to calculate pore size distributions for carbonate samples containing vugs. The results ofthis investigation indicate that, although relatively acurate results may be determined with the fastdiffusion method, the pore sizes in such samples are not all within the fast diffusion regime;therefore the slow diffusion method is more reliable for these samples. Indications of pore sizesoutside the fast diffusion regime include the stretching of the fast diffusion pore size spectra at thelargest pore sizes for samples S54, S55 and S77, and the fact that both the bulk relaxation rate andself diffusion coefficients have an obvious effect on the pore sizes calculated with the slowdiffusion method.The primary conclusion drawn from this investigation is that the inversion of NIvIR data tocalculate pore sizes is a very unstable problem. Three possible solutions are to select the type ofpore size distribution to be determined, to incorporate more constraints into the inversion, or toobtain more accurate data over the largest possible range of measurement times. Despite theinstability of the inversion, the slow diffusion results are more accurate than those of the fastdiffusion methods for vuggy carbonate samples. The results indicate the range of pore sizespresent in the sample and at least approximate the average pore size, as shown by comparing theconfocal and slow diffusion pore size spectra.The second objective of this investigation was to determine the amount of informationabout fluid distributions that can be extracted from NMR spin-lattice relaxation data on partiallysaturated samples. The shift towards lower relaxation rates and lower amplitudes as the saturationdecreases indicates two things: lower amplitudes reflect the presence of less fluid in the sample,while lower relaxation times indicate those portions of the fluid with the longest relaxation timesare being removed preferentially. The fluid with the longest relaxation times fill the largest poreswhen the pores are in the fast diffusion regime and are found in the center of the largest pores for110Carbonate (S55) Pore Size Spectra0.4Slow Diffusion SpectraConfocal Spectra0.3 Fast Diffusion Spectra0.—C.) I,3- I0.2J0 I I I0.1 I ••]•• .f: I0.0iO-7 i06 io5 io4 io3Pore Radius (m)Carbonate (S77) Pore Size Spectra0.4Slow Diffusion SpectraConfocal Spectra____I IO Fast Diffusion SpectraII Iii III I0.2 IE ii-II Iii Io ii I> I Iii II I I0.0 ..-l,.106 io lo ioPore Radius (m)Figure 3.22: Comparisons of the pore size spectra calculated from the NMR data with the fast andslow diffusion methods with the spectra determined using the confocal laser scanning microscope.111samples with pores in the slow diffusion regime. This fact, coupled with the overall decrease inrelaxation times with saturation, indicates that although a range of pore sizes are being drainedconcurrently, the fluid is removed from the largest pores more easily.The anomalous increase in Ti spectra at low saturations for samples S41, S42, S44, andS55, is probably due to the method of de-saturation used. The evaporation (drainage) of fluidfrom a sample tends to empty the largest pores first. However, even the largest pores in the centerof the sample will not drain until the fluid within the pore is in contact with air, necessitating theemptying of many pores in the outer area of the sample before pores in the center will begin toempty. Also, if the accessibility to a large pore is limited by narrow pore throats, the large porewill empty only after the pore throats. Therefore the increase in the relaxation times at lowsaturations could indicate either the trapping of fluid in large pores due to their distance from theedges of the sample or limited accessiblity due to very narrow pore throats connecting some of thelarger pores to the remainder of the pore network.The correlation of saturations with M(O) was quite accurate, and this result, along with theexpected emptying of large pores first indicates that NMR could be used to determine more detailedinformation on fluid distributions if theoretical representations of partially saturated pores weredeveloped.112CHAPTER FOURSummaryINTRODUCTIONThe primary objective of this thesis was to explore two new methods of characterizingporous media: confocal laser scanning microscopy (CLSM) and nuclear magnetic resonance(NMR). The benefits of these two methods include their non-destructive nature and the relativespeed of the measurements. The primary advantage of the confocal microscope over othermicroscopes is its ability to produce serial sections by optically sectioning a sample, rather thanphysically sectioning it. Several advantages of NMR include the possible use of wet samples, thedirect relation between pore sizes and relaxation times, and the fact that it may be possible todetermine fluid distributions from decay curves of partially saturated media.CONFOCAL MICROSCOPYThree samples, a sandstone (Berea 100) and two limestones (S55 and S77), werecharacterized using the CLSM. The results for Berea 100 compared favourably with a Berea 300pore size distributiont from the literature, with an average pore size of 25.3 p.m determined usingthe CLSM compared to 18.0 p.m determined for Berea 300 employing capillary pressure methods(Crocker et a!., 1983). The pore size spectra determined using these two methods are also verysimilar in shape, as shown in Figure 4.1. Therefore the pore size spectra determined for samplesS55 and S77 are assumed to be accurate representations of the actual pore sizes. The pore sizespectra for samples S55 and S77 are shown in Figure 4.2. The average pore sizes determined forthese two samples was 17.5 and 90.2 p.m for S55 and S77 respectively.11340—100Mean 18.0MmMode 10.0Mm.2 0.40C)Ia)E—0>1 2 4 10 20 40 100Pore diameter (jzml16 ImagesII. 200.08_0.060.040.020.00. I106 io5 ioPore Radius (m)Figure 4.1: Berea 100 pore size distributions. The upper distribution is a capillaiypressure pore size distribution (after Crocker et aL, 1983) and the lower is thedistribution determined from 16 confocal laser scanning microscope images.11420 Images_____L______Figure 4.2: Carbonate (S55 and S77) cumulative pore size distributions. The upperplot represents the summation of the distributions of 20 S55 images, the lower of13 S77 images.0.120.080.040.001060.080.060.040.02C.—C)jC—4Cio5Pore Radius (m)1013 Images0.00•106 io-5Pore Radius (m)io3115The ability of the confocal microscope to optically section samples was used to produceserial sections of the samples in the hope that reconstruction would allow three-dimensional imageanalysis. However, the samples were fragile and therefore enough serial sections to adequatelyrepresent the pores could not be collected. Partial sets of serial sections were collected, an exampleof which is shown in Figure 2.16, and the quality of these images indicate that with furtherdevelopment of the methods used to produce the pore casts, confocal laser scanning microscopycould be a powerful tool for three-dimensional characterization of porous media. The primarydifficulty is one of computing power. The images are 384 K in size, and for detailed work, 100 ormore serial sections could be collected at spacings as small as 0.2 ELm. The storage and analysis ofsuch data would require a powerful computer with a large memory to allow reconstruction andanalysis of the images.NMRThis thesis focusses on several approaches to NMR which have not previously been used:direct comparison of fast diffusion and slow diffusion results for samples with pore sizes in bothregimes, focussing on carbonates to determine the accuracy with which pores in carbonates may bedetermined, and measuring partially saturated samples to investigate the possibility of determiningfluid disthbutions with NMR.Determination of Pore SizesThe fast diffusion method has, until recently, been the only method used to calculate poresizes from NMR data. Despite evidence that pores in some samples used for other studies werenot in the fast diffusion regime, only Davies et al. (1990) have applied slow diffusion theory to thisproblem. Rather than utilizing the more complex slow diffusion relations, most investigations116either assumed or theoretically proved that the effect of the larger pores on the calculated pore sizedistributions was minimal.The pore sizes determined using the fast and slow diffusion methods are very similar tothose determined using the CLSM , as shown in the pore size spectra of S55 and S77 in Figure3.22. The inaccuracy of the fast diffusion pore size spectra lies primarily in the stretching of thedistribution at the highest pore sizes due to the fact that the equation used to calculate pore sizesfrom T1 values becomes undefmed as T1 approaches Tib. This occurs due to the term:(1 1 -i‘T1 T1b’in the equation used to calculate pore sizes from T1 values with the fast diffusion method. Thesimilarity between the T1 and fast diffusion pore size spectra for samples S31 to S49 is illustratedby the comparison plots in Figure 4.3 for samples S3 1 and S39. This correspondence between thetwo spectra indicates that these samples are within the fast diffusion regime. Further support forthe assumption that samples S31 to S49 are in the fast diffusion regime and samples S54, S55 andS77 are in the slow diffusion regime is provided by the nature of the variation in the slow diffusionpore sizes when the surface sink parameter, the self diffusion coefficient, and the bulk relaxationrate are varied.Theoretically, small variations in the bulk relaxation rate and self diffusion coefficientshould not affect the calculation of pore sizes for samples in the fast diffusion regime while similarvariations in the surface sink parameter should have no effect for samples in the slow diffusionregime. Thus, the indifference of the results of samples S3 1 through S49 to variations in the bulkrelaxation rate and the randomness of their response to variations in the self diffusion coefficientindicate these samples are within the fast diffusion regime. In comparison, the pore sizescalculated for samples S54, S55 and S77 are affected by all three parameters, indicating thepresence of pore sizes in both the fast and slow diffusion regime.117Ti and Pore Size Distributions0.08Pore Sizes Ti0.060C)1:1.4 0.04a)0.0200.00 •••••108 io 106 i05 i03 10 10W’ 100Ti (ms) and Pore Radius (m)(a)Ti and Pore Size Distributions0.08Pore Sizes Ti0.060.—C)1:1.4 0.04..%j0 0.020.00•108 io7 106 io- i010 10 10_i 10°Ti (ms) and Pore Radius (m)(b)Figure 4.3 T1 and fast diffusion pore size distributions. (a) Sample S39 (b) SampleS31118The fact that samples S54, S55 and S77 contain pore sizes in both the fast and slowdiffusion regime should ensure that the slow diffusion method produces the most accurate poresize spectra for these samples. The graphs of the confocal, slow diffusion and fast diffusion poresize spectra of these two samples in Figure 3.22 indicate that this is the case only if the peak at themaximum pore size is neglected. This peak represents 60% of the total distribution and, with itincluded, average pore sizes of 35 and 558 p.m are calculated for S77 and S55 respectively.However, the inversion of theoretical data indicated that the addition of even 1% noise to atheoretical decay curve will produce a peak at the largest pore size which represents approximately30% of the total spectra. Therefore, the final peak is probably due to noise and may be neglected.The average pore sizes calculated neglecting this peak are 19.3 and 85.3 urn, both comparing veryfavourably to the confocal averages of 17.5 and 90.2 p.m for samples S55 and S77 respectively.Fluid DistributionsThe application of NIvIR to determine fluid distributions indicated that the fluid whichevaporates from the sample most easily is that portion of the fluid with the longest relaxation time.This is shown by the decrease in the T1 values in the spectra calculated at lower saturationscompared to those calculated at higher saturations. The portion of the fluid with the longestrelaxation time is found in the largest pores when the pores are in the fast diffusion regime and inthe center of the largest pores for samples with pores in the slow diffusion regime. This fact,coupled with the overall decrease in relaxation times with saturation, indicates that although a rangeof pore sizes are being drained concurrently, the fluid is removed from the largest pores moreeasily.The anomalous increase in Ti spectra at low saturations for samples S41, S42, S44, andS55, is probably due to the method of de-saturation used. The evaporation (drainage) of fluidfrom a sample tends to empty the largest pores first. However, even the largest pores in the centerof the sample will not drain until the fluid within the pore is in contact with air, necessitating the119emptying of many pores in the outer area of the sample before pores in the center will begin toempty. Also, if the accessibility to a large pore is limited by narrow pore throats, the large porewifi empty only after the pore throats do. Therefore the increase in the relaxation times at lowsaturations could indicate either the trapping of fluid in large pores due to their distance from theedges of the sample or limited accessiblity due to very narrow pore throats connecting some of thelarger pores to the remainder of the pore network.CONCLUSIONSThe ability to accurately characterize the pore space of geophysical materials is central toour understanding of the physical properties of these materials. In this study, it has been shownthat scanning confocal laser microscopy and nuclear magnetic resonance can provide muchinformation about the sizes and shapes of pores, and the nature of the fluids present. In order tofurther expand the application of these technologies to determine pore network characteristics, theactual geometrical shape of the pores must be approximated more closely, a more stable method ofinverting NMR data in the slow diffusion regime must be developed. The collection of moreaccurate data over the largest possible measurement time range should also improve the resultsdramatically. However, if this method is ever to be applied to data collected in situ, it must bedeveloped to deal with relatively noisy data which does not extend to measurement times below 20ms (Neuman and Brown, 1982). The determination of fluid distributions could also be improvedby determining the response of partially saturated pores and then using this information in theanalysis of the data.1 20REFERENCESBacallao, R., and Steizer, E.H.K.1989: Preservation of Biological Specimens for Observation in a Confocal FluorescenceMicroscope and Operational Principles of Confocal Fluorescence Microscopy;Methods in Cell Biology, v. 31, p. 437-453.Bloch, F1946: Nuclear Induction; Physical Review, v. 70, p. 460-474.Borgia, G.C., Fantazzini, P., and E. Mesini1990: Water ‘H Spin-Lattice Relaxation as a Fingerprint of PorousMedia; MagneticResonance Imaging, v. 8, p. 435-447.Bourbie, T., and B. Zinszner1984: Saturation Methods and Attenuation Versus Saturation Relationships inFontainebleau Sandstone; Fifty-Fourth International Meeting, Society ofExploration Geophysicists, Expanded Abstracts, p. 344-350.Brown, R.J.S., and I. Fatt1956: Measurements of Fractional Wettability of Oilfield Rocks by the Nuclear MagneticRelaxation Method; Petroleum Transactions, American Institue of Mining anedMetallurgical Engineers, v. 207, p. 262-264.Brownstein, K.R., and C.E. Tarr1977: Spin-Lattice Relaxation in a System Governed by Diffusion; Journal of MagneticResonance, v. 26, p. 17-24.Brownstein, K.R., and C.E. Tarr1979: Importance of Classical Diffusion in NMR Studies of Water in Biological Cells;Physical Review A, v. 19, p. 2446-2453.Brunauer, S., Deming, L.S., Deming, W.E., and Teller, E.1938: Adsorption of Gases in Multimolecular Layers; American Chemical SocietyJournal, v. 60, p. 309-319.Callaghan, P.T., Jolley, K.W., and R.S. Humphrey1983: Diffusion of Fat and Water in Cheese as Studied by Pulsed Field Gradient NMR;Journal of Colloid and Interface Science, v. 93, p. 52 1-529.Cercone, K.R., and Pedone, V.A.1987: Fluorescence (Photoluminescence) of Carbonate Rocks: Instrumental and AnalyticalSources of Observational Error; Journal of Sedimentary Petrology, v. 57(4), p.780-782.Coates, G.R., Muller, M., and G. Henderson1991: An Investigation of a New Magnetic Resonance Imaging Log; SPWLA Thirty-Second Annual Logging Symposium, Paper D.Crocker, M.E., Donaldson, E.C., and Marchin, L.M.1983: Comparison and Analysis of Reservoir Rocks and Related Clays; Society ofPetroleum Engineers Paper 11973.121Davies, S., and K.J. Packer1990: Pore-size distributions from nuclear magnetic resonance spin-lattice relaxationmeasurements of fluid-saturated porous solids. I. Theory and simulation; Journal ofApplied Physics, v. 67(6), p. 3163-3170.Davies, S., Kalam, M.Z., Packer, K.J., and F.O. Zelaya1990: Pore-size distributions from nuclear magnetic resonance spin-lattice relaxationmeasurements of fluid-saturated porous solids. II. Applications to reservoir coresamples; Journal of Applied Physics, v. 67(6), p. 3171-3176.Domenico, S.N1976: Effect of Brine-Gas Mixture on Velocity in an Unconsolidated Sand Reservoir;Geophysics, v. 41, p. 882-894.Domenico, S.N.1977: Elastic Properties of Unconsolidated Porous Sand Reservoirs; Geophysics, v. 42,p. 1339-1368.Doyen, P.M.1988: Permeability, Conductivity and Pore Geometry of Sandstone; Journal ofGeophysical Research: v. 93, p. 7729-7740.Egger, M.D., and Petran, M.1967: New Reflected-Light Microscope for Viewing Unstained Brain and Ganglion Cells;Science, v. 157, p. 305-307.Ehrlich, R., Kennedy, S.K., Crabtree, S.J., and Cannon, R.L.1984: Petrographic Image Analysis, I. Analysis of Reservoir Pore Complexes; Journal ofSedimentary Petrology, v. 54(4), p. 1365-1378.Evamy, B.D.1969: The Precipitational Environment and Correlation of Some Calcite Cements DeducedFrom Artificial Staining; Journal of Sedimentary Petrology, v. 39(2), p. 787-821.Fukushima, E. and S.B.W. Roeder1981: Experimental Pulse NMR: A Nuts and Bolts Approach; Addison-WesleyPublishing Company, Reading, Massachusetts, 539 p.Gallegos, D.P., Munn, K., Smith, D.M., and D.L. Stermer1987: A NMR Technique for the Analysis of Pore Structure: Application to Materials withWell-Defined Pore Structure; Journal of Colloid and Interface Science, v. 119(1),p. 127-140.Gallegos, D.P., and D.M. Smith1988: A NMR Technique for the Analysis of Pore Structure: Determination of ContinuousPore Size Distributions; Journal of Colloid and Interface Science, v. 122(1), p.143- 153.Gallegos, D.P., Smith, D.M., and C.J. Brinker1988: A NMR Technique for the Analysis of Pore Structure: Application to Mesoporesand Micropores; Journal of Colloid and Interface Science, v. 124(1), p. 186-198.1 22Gies, R.M.1987: An Improved Method for Viewing Micropore Systems in Rocks with the PolarizingMicroscope; Society of Petroleum Engineers Formation Evaluation, v. 2(2), p.209-214.Grayson, J.F.1956: The Conversion of Calcite to Fluorite; Micropaleontology, v. 2, p. 7 1-78.Hamilton, D.K., and Wilson, T.1982: Surface Profile Measurements using the Confocal Microscope; Journal of AppliedPhysics, v. 7, p. 5320-5322.Hedberg, S.A.1990: The NMR Response of Sandstones and its application to the detection ofHydrocarbon Contamination; B.A.Sc. Thesis, University of British Columbia,Vancouver, British Columbia, Canada.Jackson, J.A.1984: Nuclear Magnetic Resonance Well Logging; The Log Analyst, September-October,p. 16-30.Kenyon, W.E., Day, P.1., Straley, C., and J.F. Willemsen1988: A Three-Part Study of NMR Longitudinal Relaxation Properties of Water-SaturatedSandstones; SPE Formation Evaluation, p. 622-633.Kenyon, W.E., Howard, J.J., Sezginer, A., Straley, C., Matteson, A., Horkowitz, K., andR.Ehrlich1989: Pore-size Distribution and NMR in Microporous Cherty Sandstones; Society ofProfessional Well Log Analysts Thirtieth Annual Logging Symposium, Paper LL.Kleinberg, R.L., and Horsfield, M.A.1990: Transverse Relaxation Processes in Porous Sedimentary Rock; Journal of MagneticResonance, v. 88, p. 9-19.Knight, R.1991: Hysteresis in the Electrical Resistivity of Partially Saturated Sandstones;Geophysics, v. 56(12), p. 2139-2147.Knight, R., and R. Nolen-Hoeksema1990: A Laboratory Study of the Dependence of Elastic Wave Velocities on Pore ScaleFluid Distribution; Geophysical Research Letters, v. 17(10), p. 1529-1532.Knight, R., and A. Nur1987: Geometrical Effects in the Dielectric Response of Partially Saturated Sandstones;The Log Analyst, v. 28, p. 513-5 19.Knispel, R.R.1981: Review of Two Models for Water Proton Relaxation in Tissue; Bulletin of MagneticResonance, v. 3, p. 104-108.Lin, C., and Hamasaki, J.1983: Pore Geometry: A New System for Quantitative Analysis and 3-D display; Journalof Sedimentary Petrology, v.53, p. 670-672.1 23Longeron, D.G., Arguad, M.J., and J.P. Feraud1989: Effect of Overburden Pressure and the Nature and Microscopic Distribution ofFluids on Electrical Properties of Rock Samples; Society of Petroleum EngineersFormation Evaluation, v. 4, p. 194-202.Loren, J.D., and J.D. Robinson1970: Relations Between Pore Size Fluid and Matrix Properties, and NIvIL Measurements;Society of Petroleum Engineers Journal, v. 249, p. 268-278.Lukosz, W.1966: Optical Systems with Resolving Powers Exceeding the Classical Limit; Journal ofthe Optical Society of America, v. 56, p. 1463-1472.Marshall, G.E., Konstas, A.G., and Lee, W.R.1990: Ultrastructural distribution of collagen yypes 1-VI in aging human retinal vessels;British Journal of Opthamology, v. 74(4), p. 228-232.Masters, B.R., Kriete, A., and Kukulies, J.1991: Fluorescence Confocal Microscopy of the Living Eye; Scanning, v. 13, p. 1-16.Menon, R.E., MacKay, A.L., Hailey, J.R.T., Bloom, M., Burgess, A.E., and J.S. Swanson1987: An NMR determination of the physiological water distribution in wood duringdrying; Journal of Applied Polymer Science, v. 33, p. 1141-1155.Minsky, M.1961: Microscopy Apparatus; US Patent 3013467, Dec. 19 (Filed Nov. 7, 1957).Munn, K., and D.M. Smith1987: A NMR Technique for the Analysis of Pore Structure: Numerical Inversion ofRelaxation Measurements; Journal of Colloid and Interface Science, v. 119(1), p.117- 126.Neuman, C.H., and R.J.S. Brown1982: Applications of Nuclear Magnetism Logging to Formation Evaluation; Journal ofPetroleum Technology, v. 34(12), p. 2853-2862.Nutting, P.G.1929: Some physical problems in oil recovery; Oil and Gas Journal, v. 28, p. 44-45.Parker, K.J. and C. Rees1972: Pulsed NMR Studies of Restricted Diffusion I. Droplet Size Distributions inEmulsions; Journal of Colloid and Interface Science, v. 40(2), p. 206-2 18.Penner, D.G.1957: Facies and porosity relationships in the Mississippian Elkton carbonate cycle ofsouthwestern Alberta; Alberta Society of Petroleum Geologists Journal, v. 5(5), p.101-104.Pierce, C.1953: Computation of Pore Sizes from Physical Adsorption Data; Journal of PhysicalChemistry, v. 57, p. 149-152.1 24Pittman, E.D., and Duschatko, R.W.1970: Use of pore casts and scanning electron microscopy to study pore geometry;Journal of Sedimentary Petrology, v. 40, p. 1153-1157.Rassineux, F., Beaufort, D., Meunier, A., and Bouchet, A.1987: A Method of Coloration by Fluorescein Aqueous Solution for Thin-SectionMicroscopic Observation; Journal of Sedimentary Petrology, v. 57, p. 782-783.Robinson, J.D., Loren, J.D., Vajnar, E.A., and D.E. Hartman1974: Determining Residual Oil with the Nuclear Magnetism Log; Journal of PetroleumTechnology, v. 257, p. 226-236.Russ, J.C.1990: Computer Assisted Microscopy: The Measurement and Analysis of Images, PlenumPress, 453 p.Russ, J.L., Palmour III, H., and Hare, T.M.1989: Direct 3-D pore location measurement in alumina Journal of Microscopy, v. 155,p. RP1-RP2.Ruzyla, K., and Jezek, D.I.1987: Staining Method for Recognition of Pore Space in Thin and Polished Sections;Journal of Sedimentary Petrology, v. 57, p. 777-778.Schmidt, E.J., Velasco, K.K., and A.M. Nur1986: Quantifying solid-fluid interfacial phenomena in porous rocks with proton nuclearmagnetic resonance; Journal of Applied Physics, v. 59(8), p. 2788-2797.Shuman, H., Murray, J.M., and DiLullo, C.1989: Confocal Microscopy: An Overview; Biotechniques, v. 7, p. 154-163.Smart, P., and Tovey, N.K.1982: Electron microscopy of soils and sediments: techniques; Clarendon Press, Oxford,p. 61-80.Soeder, D.J.1990: Applications of Fluorescence Microscopy to Study of Pores in Tight Rocks;American Association of Petroleum Geologists Bulletin, v. 74, p. 30-40.Spurr, A.R.1969: A Low-Viscosity Epoxy Resin Embedding Medium for Electron Microscopy;Journal of Ultrastructure Research, v. 26, p. 3 1-43.Straley, C., and Minnis, M.M.1983: Epoxy Rock Replicas for Microtoming; Journal of Sedimentary Petrology, v. 53,p. 667-669.Straley, C., Morriss, C.E., Kenyon, W.E., and J.J. Howard1991: NMR in Partially Saturated Rocks: Laboratory Insights on Free Fluid Index andComparison with Borehole Logs; Society of Petroleum Well Log Analysts ThirtySecond Annual Logging Symposium, Paper CC.Tercier, P.1992: Electrical Resistivity of Partially Saturated Sandstones; M.Sc. Thesis, University ofBritish Columbia, Vancouver, British Columbia, Canada.1 25Thomas, G.E. and R.P. Glaister1960: Facies and Porosity Relationships in some Mississippian Carbonate Cycles of theWestern Canada Basin; Bulletin of the American Association of PetroleumGeologists, v. 44(5), p. 569-588.Timur, A.1969: Pulsed Nuclear Magnetic Resonance Studies of Porosity, Movable Fluid, andPermeability of Sandstones; Journal of Petroleum Technology, v. 246, P. 775-786.Vinegar, H.J., Tutunjian, P.N., Crabtree, P.T., Taffaldi, R.J., DiFoggio, R., and W.A. Edeistein1989: NMR Spectroscopy of Tight, Gypsum-Bearing Carbonates; Society of ProfessionalWell Log Analysts Thirtieth Annual Logging Symposium, Paper HH.Wardlaw, N.C.1976: Pore Geometry of Carbonate Rocks as Revealed by Pore Casts and CapillaryPressure; American Association of Petroleum Geology Bulletin, v. 60(2), p. 245-257.Whittall, K.P. and A.L. MacKay1989: Quantitative Interpretation of NMR Relaxation Data; Journal of MagneticResonance, v. 84, p. 134-152.Whittall, K.P.1991: Recovering Compartment Sizes from NMR Relaxation Data; Journal of MagneticResonance, v. 94, p. 486-492.Wicksell, S.D.1926: The corpuscle problem. Second memoir. Case of ellipsoidal corpuscles;Biometrika, v. 18, p. 152-172.Wilson, T. (ed.)1990: Confocal Microscopy, Academic Press, Great Britain, 426 p.Winsauer, W.O., Shearin, H.M. Jr., Masson, P.H., and Williams, M.1952: Resistivity of brine-saturated sands in relation to pore geometry; AmericanAssociation of Petroleum Geology Bulletin, v. 36, p. 253-277.Young, F.G.1967: Elkton Reservoir of Edson Gas Field, Alberta; Bulletin of Canadian PetroleumGeology, v. 15(1), p. 50-64.Young J.Z., and Roberts, F.1951: A Flying-Spot Microscope; Nature, v. 167, p. 231.Zimmerman, J.R. and W.E. Brittin1957: Nuclear Magnetic Resonance Studies in Multiple Phase Systems: Lifetime of aWater Molecule in an Adsorbing Phase on Silica Gel;Journal of PhysicalChemistry, v. 61(6), p. 1328-1333.1 26APPENDIX 1BioRad MRC-500Image File FormatThese picture files are standard MS-DOS files which consist of a 76 byte header, a pixelarray of dimensions NX x NY x Npic, and a series of note records of 96 bytes each.Header StructureContents0-1 NX - X pixel dimension of the image2-3 NY - Y pixel dimension of the image4-5 Npic - number of image frames in the file6-7 Contrast ramp 0 (pixel value for full black)8-9 Contrast ramp 255 (pixel value for peak white)10-13 Ignore14-15 1 if byte pixels, 0 if word pixels16-49 Ignore50-51 1 if merged image, 0 if unmerged image52-53 Low byte used for selecting RGB colours (7=white)54-55 Set to 12345, used to check for valid picture file56-57 Contrast ramp 0 for second picture of merged pair58-59 Contrast ramp 255 for second picture of merged pair60-6 1 Low byte used for RGB colour of second picture of merged pair64-65 lens power66-69 Magnification factor, real variable - see note below70-75 IgnoreNOTE: the length of a line of pixels in microns is computed as:((# pixels)x(Scale factor))/((lens mag.)x(mag. factor))where the scale factor is the pixel to micron ratio which is set up in the setup menu, and themagnification factor is the value stored in the file header.Pixel Array:The pixel array contains Npic images with X and Y pixel dimensions of NX and NY. Thisarray is either NX x NY x Npic bytes in size or, if word (16 bit) pixels rather than the usual byte (8bit) pixels are used, the array is 2 x (NX x NY x Npic) bytes in size.1 27-L I\)I.—JcI9?-‘.0‘IiI-C—to;10 I0IICD-.I—-CDCD0—‘.+fCD-4•E‘. cJ-—. ic. ‘-‘D• 00 -f- CD E Ii r.crQ O 0 ri CDAPPENDIX 2Image AnalysisProgramsBINREP.FCC FEBUARY, 1992 ALICE CHAPMANC ROCK PHYSICS LAB, UBCCC THIS PROGRAM WILL READ IN CONFOCAL IMAGE FILES WRITTENC IN THE FORMAT 10013 AND AN ARRAY SIZE OF 768 BY 512C PIXELS. THE INTENSITY DISTRIBUTIONS OF THE IMAGES ARE DEC TERMINED AND WR11TEN TO AN OUTPUT FILE, A PORE/SOLIDC DIVISION IS DONE TO PRODUCE BINARY IMAGES, AND THE SMALLC HOLES IN THE PORES ARE CLOSED USING A DILATION/EROSIONC SEQUENCE. FINALLY, SMALL PORES WHICH ARE ASSUMED TO BEC ARTIFACTS OF THE CHOSEN PORE/SOLD CUTOFF VALUE AREC REMOVED FROM THE IMAGES USING AN EROSION/DILATION SEQUENCE.CC MODIFIED FROM: BIN.FCC***************************CC THE VARIABLES ARE DECLARED.CINTEGER*2 IMG(530,780)INTEGER POR,SOL,INTNS(256),DIL,ERO,NREPCHARACTER*20 INP(20),OUTP(20)CC***************************CC VARIABLE DESCRIPTION:CC IMG(530,780): THE CONFOCAL IMAGE IS PLACED IN THIS ARRAYC WITH THREE ROWS OF SOLIDS SURROUNDING ITC POR: THE NUMBER WHICH PIXELS WITHIN THE PORES OF THE IMAGEC WILL BE SET TO (255)C SOL: THE NUMBER WHICH PIXELS WiTHIN THE SOLID PORTION OF THEC 11VIAGE WILL BE SET TO (0)C INTNS(256): A RECORD OF THE NUMBER OF PIXELS IN THE ORIGINALC IMAGE WHICH HAD A FLUORESCENT INTENSiTY OF 0 TO 255C DIL: THE NUMBER WHICH SOLID PIXELS IN CONTACT WITH THE PORESC WILL BE SET TO DURING THE DILATION SEQUENCEC ERO: THE NUMBER WHICH PORE PIXELS IN CONTACT WITH THE SOLD)C WILL BE SET TO DURING THE EROSION SEQUENCEC NREP: THE NUMBER OF IMAGES WHICH WILL BE ALTERED WITH THISC PROGRAMC INP: THE NAMES OF THE INPUT FILES (USUALLY THE NAMES OF THEC MODELS)C OUTP: THE NAMES OF THE OUTPUT FILESCC***************************C1 29C ThE WIDTH (NX) AND LENGTH (NY) OF THE CONFOCAL IMAGE ARRAYC AND THE WORKING ARRAY (NTX) AND (NTY) ARE SET.CNX = 512NY = 768NTX = NX + 6NTY = NY +6CC THE INTENSITY OF FLUORESCENT LIGHT RESPONSE WHICH CORRESPONDSC TO THE CUT OFF BETWEEN PORE AND SOLID IN THE ORIGINAL IMAGE ISC SET BY THE USER.CWpJTE(9,*) “WHAT CUT OFF WOULD YOU LIKE FOR THE SOLUYr’READ(9,*) NSOLCC THE NUMBER OF IMAGES WHICH WILL BE ALTERED IS INiTIALIZED TO ZERO,C THEN THE USER IS REQUESTED TO ENTER THE DESIRED NUMBER.C11 NREP=0WPJTE(9,*)”PLEASE ENTER THE NUMBER OF IMAGES& YOU WISH TO HAVE ALTERED (0 <INT < 100)”READ(9,*)NREPCC THE NUMBER OF IMAGES TO BE ALTERED IS CHECKED TO ENSURE THATC iT IS LESS THAN 100 AND IS NOT NEGATIVE. IF IT DOES NOT MEETC THESE REQUIREMENTS, A MESSAGE IS WRITTEN FOR THE USER, ANDC THEY ARE ASKED TO ENTER AN ACCEPTABLE VALUE.CIF ((NREP.LT.1).OR.(NREP.GT.99)) THENIF (NREP.LT.1) THENwPJTE(9,*)”THIs NUMBER IS NEGATIVE!”ELSEWRITE(9,*)”THIS NUMBER IS TOO LARGE!”ENDIFGOTO 11ENDIFCC THE NAME OF THE INPUT AND OUTPUT FILES ARE REQUESTED FROM THEC USER.CDO 5 I=1,NREPwR1TE(9,*)”EITER THE IMAGE FILE NAME”READ(9,*)INP(I)WRITE(9,*) “ENTER THE OUTPUT FILE NAME”READ(9,*)OUTP(I)5 CONTINUECC THE LOOP WHICH WILL ALTER AND WRiTE OUT EACH IMAGE ISC BEGUN.CDO 200 NINC = 1,NREPCC THE NINC (FIRST,SECOND...) INPUT AND OUTPUT FILES ARE OPENED.1 30COPEN(5 ,FILE=INP(NINC))OPEN(6,FILE=OUTP(NINC))CC THE TWO ARRAYS IN THE PROGRAM ARE IMTIALIZED TO ZERO.CDO 10I=1,NTXDO 10 J=1,NTYIMG(I,J) = 010 CONTINUEDO 20 1=1,2561NTNS(I) =020 CONTINUECC THE IMAGE IS READ IN FROM THE IMAGE FILE.CREAD(5,91 )NTEST,NTEST291 FORM.AT(/////////18/18)DO 30 I=4,NX÷3WRITE(9,*)”I = 11,1READ(5,92)(IMG(I,J),J=4,NY÷3)30 CONTINUE92 FORMAT(10013)CC THE INTENSITY OF EACH PIXEL IN THE ORIGINAL IMAGE ISC RECORDED IN THE ARRAY INTNS AND WRHTEN TO THE FILEC “INTNS”INP WHICH WAS OPENED FOR THAT PURPOSE.COUTP(NINC) = “INTNS”/IINP(NINC)OPEN(7,FILE = OUTP(NTNC))CDO 40 I=4,NX+3DO 40 J=4,NY+3ITMP = IMG(I,J) + 1INTNS(ITMP) = INTNS(1TMP) + 140 CONTINUECWRITE(7,93)93 FORMAT(”INTENSITY DISTRIBUTION”)CDO 501=1,256WRITE(7,94)I,INTNS(I)50 CONTINUE94 FORMAT(2110)CCLOSE(7)CC THE PORE/SOLID DIVISION IS PERFORMED TO PRODUCE A BINARYC IMAGE OF ZERO’S (SOLID) AND 255’S (PORE).CSOL=0POR = 255C131DO 60 1=4,NX+3DO 60 J=4,NY+3IF (IMG(I,J).GT.NSOL) THENIMG(I,J) = PORELSEIMG(I,J) = SOLENDIF60 CONTINUECC NOW AN EROSION/DILATION SEQUENCE IS PERFORMED TO REMOVE ANYC SMALL PORES WHICH ARE SIMPLY ARTIFACTS OF THE PROCESS OFC CHANGING THE IMAGE TO A BINARY IMAGE.CERO = 245CC ANY PORE PIXELS IN CONTACT WITH A SOLID ARE RE-NUMBERED TOC THE ERO VALUE OF 245.CDO 70 I=4,NX+3DO 70 J=4,NY÷3IF (IMG(I,J).EQ.SOL) THENIF (IMG(I÷1,J).EQ.POR) IMG(I-i-1,J) = EROIF (IMG(I-1,J).EQ.POR) IMG(I-1,J) = EROIF (IMG(I,J+1).EQ.POR) IMG(I,J-i-1) = EROIF (IMG(I,J-1).EQ.POR) IMG(I,J-1) = ERO14 ENDIF70 CONTINUECC PIXELS WHICH WERE CHANGED TO ERO AN]) ARE STILL IN CONTACTC WITH THE PORE ARE RE-SET TO POR(255). EACH PIXEL IN THEC IMAGE IS CHECKED AND THE IMAGE IS GONE THROUGH ONLY ONCE TOC ENSURE THAT SMALL ROUGHNESSES ON THE PORES EXTENDING INTOC THE SOLID ARE SMOOTHED.CDO 80 I=4,NX+3DO 80 J=4,NY+3IF (IMG(I,J).EQ.POR) THENIF (IMG(I+1,J).EQ.ERO) IMG(I+1,J) = PORIF (IMG(I-1,J).EQ.ERO) IMG(I-1,J) = PORIF (IMG(I,J+1).EQ.ERO) IMG(I,J+1) = PORIF (IMG(I,J-1).EQ.ERO) IMG(I,J-1) = POR15 ENDIF80 CONTINUECC ANY ALTERED PIXELS WHICH ARE EITHER COMPLETELY SURROUNDED BYC SOLID PIXELS OR ARE PART OF EXTENWE ROUGHNESS ON THE POREC SURFACE ARE STILL SET TO ERO AN]) ARE NOW IRREVOCABLYC RE-SET TO SOLID PIXELS (0).C16 DO 90 I=4,NX÷3DO 90 J=4,NY+3IF (IMG(I,J).EQ.ERO) THENIMG(I,J) = SOLENDIF1 3290 CONTINUECC NOW A DILATION/EROSION SEQUENCE IS PERFORMED TO CLOSE ANYC SMALL HOLES IN THE PORES.CDIL = 250CC ANY SOLID PIXELS IN CONTACT WITH A PORE ARE RE-NUMBERED TOC THE DIL VALUE OF 250.CDO 100 I=4,NX+3DO 100 J=4,NY÷3IF (IMG(I,J).EQ.POR) THENIF (IMG(I+1,J).EQ.SOL) IMG(I+1,J) = DILIF (IMG(I-1,J).EQ.SOL) IMG(I-1,J) = DILIF (IMG(I,J+1).EQ.SOL) IMG(I,J+1) = DILIF (IMG(I,J-1).EQ.SOL) IMG(I,J-1) = DILENDIF100 CONTINUECC PIXELS WHICH WERE CHANGED TO DIL AND ARE STILL IN CONTACTC WITH THE SOLID ARE RE-SET TO SOL(0). EACH PIXEL IN THEC IMAGE IS CHECKED AND THE IMAGE IS GONE THROUGH ONLY ONCE,C THUS FILLING SMALL CHANNELS INTO THE PORE SPACE AN]) GENERALLYC ‘SMOOTHING’ THE BINARY IMAGE.CDO 110 I=2,NX+5DO 110 J=2,NY+5IF (IMG(I,J).EQ.SOL) THENIF (IMG(I+1,J).EQ.DIL) IMG(I+1,J) = SOLIF (IMG(I-1,J).EQ.DIL) IMG(I-1,J) = SOLIF (IMG(I,J+1).EQ.DIL) IMG(I,J+1) = SOLIF (IMG(I,J-1).EQ.DIL) IMG(I,J-1) = SOL12 ENDIF110 CONTINUECC ANY ALTERED PIXELS WHICH ARE COMPLETELY SURROUNDED BY POREC PIXELS ARE STILL SET TO DIL AND ARE NOW IRREVOCABLY RE-SETC TO PORE PIXELS (255).C13 DO 120 I=4,NX+3DO 120 J=4,NY+3IF (IMG(I,J).EQ.D1L) THENIMG(I,J) = PORENDIF120 CONTINUECC THE MODIFIED IMAGE IS NOW WREN INTO THE OUTPUT FILE.CDO 130 I=4,NX--3WRITE(6,95)(IMG(I,J),J=4,NY+3)130 CONTINUE95 FORMAT(10014)C1 33C THE INPUT AND OUTPUT FILES ARE CLOSED.CCLOSE(5)CLOSE(6)CC THE NEXT STAThMENT IS THE END OF THE LOOP WHICHC IS REPEAThD FOR EACH IMAGE.C200 CONTINUECSTOPEND1 34IMAGANAL.FCC MARCH, 1992 ALICE CHAPMANC ROCK PHYSICS LAB, UBCCC THIS PROGRAM READS iN BINARY CONFOCAL IMAGE FILESC PRODUCED BY THE PROGRAM BINIMGR.F. THESE IMAGES AREC MANIPULATED TO ALTER ALL PIXELS WiTHIN PORE ThROATSC TO VALUES OF SEPR (20). A WATERSHED ALGORiTHM ISC THEN USED TO DETERMINE THE POINTS WITHIN THE PORESC AND PORE THROATS WHICH ARE FARTHEST FROM THE EDGE OFC THE PORE OR PORE THROATS. THESE PORE AND PORE THROATC CENTERS ARE THEN USED TO DETERMINE THE SIZE OF THEC PORES AND PORE THROATS. FINALLY, THE OVERALL SURFACEC AREA AND VOLUME OF THE PORES iN THE IMAGE IS DETERMINEDC AN]) THE POROSITY (RATIO OF THE VOLUME OF PORES TOC TOTAL VOLUME OF THE IMAGE) AND SURFACE AREA TO VOLUMEC RATIO OF THE PORES IS CALCULATED.CCC***************************CC THE VARIABLES ARE DECLARED.CINTEGER*2 IMG(530,780)INTEGER POR,SOL,SEPR,BDRY,PRDIST( 101)INTEGER NPR,NTHR,IMAX,THDIST( 101),MXRADREAL POROS,PSV,TSV,THSVCHARACTER*20 INP(20),OUTP(20),MDL(20),TMPCC***************************CC VARIABLE DESCRIPTION:CC IMG(530,780): THE CONFOCAL IMAGE IS PLACED IN THIS ARRAYC WITH THREE ROWS OF SOLIDS SURROUNDING ITC POR: THE NUMBER WHICH PIXELS WITHIN THE PORES OF THEC IMAGE WILL BE SET TO (255)C SOL: THE NUMBER WHICH PIXELS WITHIN THE SOLD) PORTION OFC THEIMAGEWILLBESETTO(0)C SEPR: THE NUMBER WHICH PIXELS IN THE SMALLEST PART OF THEC PORE THROATS WILL BE SET TO TO SEPARATE THE PORESC BDRY: THE INITIAL NUMBER WHICH PORE PIXELS BORDERING THE SOLD)C ARE SET TO IN THE WATERSHED ALGORiTHM. PORE PIXELS WHICHC BORDER THOSE WHICH ARE SET TO BDRY ARE THEN SET TO BDRY+1C ETC, ETC.C CrR: THE NUMBER TO WHICH THE PIXEL AT THE CENTRE OF A POREC OR PORE THROAT IS SETC NREP: THE NUMBER OF IMAGES WHICH WILL BE ALTERED WITH THISC PROGRAMC PRDIST: A DISTRIBUTION OF THE PORE SIZES IN THE IMAGE. PRDIST(I)C WILL BE A NUMBER (J) INDICATING J PORES OF RADI1JS 1-0.5C WERE PRESENT IN THE IMAGEC THDIST: A DISTRIBUTION OF THE PORE THROAT SIZES IN THE IMAGE.1 35C THDIST(I) WILL BE A NUMBER (J) INDICATING J PORE THROATSC OF RADIUS 1-0.5 WERE PRESENT IN THE IMAGEC NPR: THE TOTAL NUMBER OF PORES WHOSE SIZE WAS MEASUREDC NTHR:THE TOTAL NUMBER OF PORE THROATS WHOSE SIZE WAS MEASUREDC MXRAD:THE INTEGER REPRESENTING THE SIZE OF THE LARGEST PORE SIZEC MEASURED. MXRAD = RADIIJS+0.5C POROS: THE POROSITY OF THE IMAGE: IE THE RATIO OF THE NUMBERC OF PORE AND PORE THROAT PIXELS IN THE IMAGE TO THE TOTALC NUMBER OF PIXELS IN THE IMAGE.C PSV: THE SURFACE TO VOLUME RATIO (PSA/PVOL) OF THE PORESC THSV: THE SURFACE TO VOLUME RATIO (THSA1HVOL) OF THE POREC THROATSC TSV:THE OVERALL SURFACE TO VOLUME RATIO OF THE VOII) SPACE IN THEC IMAGE: (PSA + THSA)/(PVOL + THVOL)C INP: THE NAMES OF THE INPUT FILESC OUTP: THE NAMES OF THE OUTPUT FILESC MDL: THE NAME OF THE MODELS, WHICH ARE USED TO NAME SOMEC OUTPUT FILESCC***************************CC THE WIDTH (NX) AND LENGTH (NY) OF THE CONFOCAL IMAGE ARRAYC AND THE WORKING ARRAY (NTX) AND (NTY) ARE SET.CNX = 512NY = 768NTX = NX + 6NTY = NY +6CC THE NUMBER OF IMAGES WHICH WILL BE ALTERED IS INITIALIZED TOC ZERO, THEN THE USER IS REQUESTED TO ENTER THE DESIRED NUMBER.C15 NREP=0WRITE(9,*)”PLEASE ENTER THE NUMBER OF IMAGES& YOU WISH TO HAVE ALTERED (0 <INT < 100)”READ(9,*)NREPCC THE VALUE OF NREP IS CHECKED TO ENSURE THAT IT IS LESS THAN 100C AND IS NOT NEGATIVE. IF IT DOES NOT MEET THESE REQUIREMENTS,C A MESSAGE IS WRITTEN FOR THE USER, AND THEY ARE ASKED TOC ENTER A REASONABLE VALUE.CIF ((NREP.LT.1).OR.(NREP.GT.99)) THENIF (NREP.LT.1) THENWRITE(9,*)”TBIS NUMBER IS NEGATIVE!”ELSEWRITE(9,*)”THIS NUMBER IS TOO LARGE!”ENDIFGOTO 15ENDIFCC THE NAME OF THE MODEL AND THE INPUT AND OUTPUT FILES AREC REQUESTED FROM THE USER.C1 36DOS I=1,NREPwPJTE(9,*)”ENTER THE IMAGE FILE NAME”READ(9,*)INP(I)WPJTE(9*)”ENTER THE OUTPUT FILE NAME”pEAD(9,*)OUTp(I)WRITE(9,*)”ENTER THE NAME OF THE MODEL”PEAD(9,*)MDL(I)5 CONTINUECC THE LOOP WHICH WILL ALTER AND WRITE OUT EACH IMAGE ISC BEGUN.CDO 380 NINC = 1,NREPCC THE NINC (FIRST,SECOND...) INPUT FILE IS OPENED.COPEN(5 ,FILE=INP(NINC))CC THE ARRAYS IN THE PROGRAM ARE IN1TIALIZED TO ZERO.CDO 10 I=1,NTXDO 1OJ=1,NTYIMG(I,J) = 010 CONTINUECC THE IMAGE IS NOW READ IN FROM THE INPUT FILE.CWRITE(9,*)”IMAGE “,NINC,” IS BEING ANALYSED”DO 30 I=4,NX+3READ(5,9 1)(IMG(I,J),J=4,NY+3)30 CONTINUE91 FORMAT(100I4)CC THE PORE/SOLID, SEPARATION, CTh, AND BOUNDARY VALUES FOR THEC IMAGE FILE ARE INITIALIZED.CPOR = 255SOL=0BDRY =50SEPR=30IMAX =63CTR=200CC THE “WATERSHED ALGORITHM IS NOW APPLIED. THIS SETS THE VALUESC OF THE PIXELS WITHIN THE PORES TO PROGRESSIVELY LOWER OR HIGHERC VALUES DEPENDING ON THEIR PROXIMITY TO THE SOLID.CDO 40 I=50,TMAXITMP = IDO 50 J= 3, NX+4DO 50 K= 3, NY÷41 37IF (I.EQ.50) THENIF (IMG(J,K).LT.SEPR) THENIF (IMG(J+i,K).EQ.POR) IMG(J÷1,K) = ITMPIF (IMG(J-i,K).EQ.POR) IMG(J-1,K) = ITMPIF (IMG(J,K÷1).EQ.POR) IMG(J,K+i) = ITMPIF (IMG(J,K-1).EQ.POR) IMG(J,K-1) = ITMPENDIFELSEIF (IMG(J,K).EQ.(ITMP-1)) THENIF (IMG(J+l,K).EQ.POR) IMG(J÷l,K) = ITMPIF (IMG(J-1,K).EQ.POR) IMG(J-1,K) = ITMPIF (IMG(J,K+1).EQ.POR) IMG(J,K-i-1) = ITMPIF (IMG(J,K-1).EQ.POR) IMG(J,K-1) = ITMPENDIFENDIF50 CONTINUE40 CONTINUEDO 360 MAX = 51,IMAX,4SEPR =30CC ALL PIXELS SET TO VALUES BETWEEN 50 AND MAX-i ARE RE-SETC TO SEPR.CDO 130 I=3,NX+4DO 130 J=3,NY+4IF (IMG(I,J).GT.SEPR) THENIF (IMG(I,J).LT.MAX) THENIMG(I,J) = SEPRENDIFENDIF130 CONTINUECC NOW THE WATERSHED ALGORITHM IS RE-APPLIED TO BUILD THE PORESC UP BASED ON THE REMAINING LAYER OF MAX’S IN THE PORES.CNTMP =48DO 140 I=MAX-1,NTMP,-1NCHK=I+i1TMP IDO 150J=3,NX+4DO 150 K= 3, NY-1-4IF (IMG(J,K).EQ.NCHK) THENIF (IMG(J+1,K).EQ.SEPR) IMG(J+1,K) = ITMPiF (IMG(J-1,K).EQ.SEPR) IMG(J-1,K) = ITMPIF (IMG(J,K+1).EQ.SEPR) IMG(J,K+1) = ITMPIF (IMG(J,K-1).EQ.SEPR) IMG(J,K-1) = ITMPENDIF150 CONTINUE140 CONTINUECC THE PIXELS NEWLY ALTERED TO THE PREVIOUS VALUE OF SEPR (30)C ARE CHECKED TO SEE WHAT THEIR SURROUNDINGS ARE. IF THERE IS1 38C APOREONONESIDEANDEITHERAPORETHROATORSOLIDONTHEC OTHER, THE PIXEL IS SET TO 15. IF ThERE WAS NOT A PORE THROATC OR SOLID PRESENT, iT IS SET TO THE CURRENT SEPR VALUE, 20.C ALSO, IF THERE IS ONLY EiTHER SOLID OR PORE THROAT CLOSE BY,C THE PIXEL IS SET TO 20.CSEPR =20DO 160 I=3,NX-i-4DO 160 J=3,NY+4IF (IMG(I,J).EQ.30) THENNSOL=0NPOR=0Nil =0DO 170 K=I-3,I+3IF (IMG(K,J).GT.47) THENNPOR= 1ELSEIF (IMG(K,J).LT.1 1) THENNSOL=iELSEIF(IMG(K,J).EQ.11) Nil = 1ENDIFENDIF170 CONTINUEIF (NPOR.EQ.1) THENIF (NSOL.NE.1).OR.(Nll.EQ.1) THENIMG(I,J) = SEPRGOTO 160ELSEIMG(I,J) = 15ENDIFELSEIF (NSOL.EQ.1).OR.(Nll.EQ.1) THENIMG(I,J) = SEPRENDIFENDIFNPOR=0DO 180 K=J-3,J+3IF (IMG(I,K).GT.47) THENNPOR=lELSEIF (IMG(I,K).LT.ll) THENNSOL=NSOL + 1ELSEIF(IMG(K,J).EQ.1l) Nil = Ni 1+1ENDIFENDIF180 CONTINUEIF (NPOR.EQ.1) THENIF (NSOL.NE.i).OR.(Ni1.NE.l) THENIMG(I,J) = SEPRELSEIMG(I,J) = 15ENDIF1 39ELSEIF (NSOL.GT.1) THENIMG(I,J) = 15ENDIFIF (N11.GT.1) THENIMG(I,J) = SEPRENDIFENDIF12 ENDIF160 CONTINUECC PIXELS SET TO 30 ARE NOW CHECKED SO THAT ANY IN CONTACTC WITH PiXELS SET TO SEPR ARE ALSO SET TO SEPR. AS THEC SURROUNDING OF EACH SEPR PIXEL IS CHECKED, THE PIXEL ISC SET TO 25 SO THAT iT WILL ONLY BE CHECKED ONCE.CDO 190 L=1,200NSEP0DO 200 I=3,NX+4DO 200 J=3,NY-i-4IF (IMG(I,J).EQ.SEPR) THENIMG(I,J) = 25DO 210 K=I-1,I-i-1,2IF (IMG(K,J).EQ.30).OR.(IMG(K,J).EQ.15) THENIMG(K,J) = SEPRNSEP = 1ENDIF210 CONTINUEDO 220 K=J-1,J+1,2IF (IMG(I,K).EQ.30).OR.(IMG(I,K).EQ. 15) THENIMG(I,K) = SEPRNSEP = 1ENDIF220 CONTINUEENDIF200 CONTINUEIF (NSEP.EQ.0) GOTO 185190 CONTINUE185 CONTINUECC PIXELS CHANGED TO 25 TO SPEED THE PREVIOUS LOOP AREC RE-SET TO 20.CDO 195 I=4,NX÷3DO 195 J=4,NY+3IF (IMG(I,J).EQ.25) IMG(I,J)=20195 CONTINUECC AFTER EACH LOOP IN WHICH MAX CHANGES, THE PIXELS WHICH WEREC PREVIOUSLY SET TO BE PORE THROATS ARE RE-SET TO 11. IF MAXC IS GREATER THAN 60, ThEN ANY NEWLY CREATED PORE-THROATS INC CONTACT WITH AN OLD PORE THROAT IS RE-SET TO 10, IE: THEYC ARE NOT ALLOWED TO REMAIN PORE THROATS. ALSO, ANY PIXELSC IN CONTACT WITH PIXELS SET TO 10 ALREADY ARE SET TO 10.140CIF (MAX.GT.60) THENM1=2M2= 1INC=-1ELSEM1=1M2=2INC= 1ENDIFDO 485 M=M1,M2,INCrRfl’E(9,*)”M =N=9+MN3 = NN2= 10+2*MIF (MAX.GT.60) THENN3=10ENDIFDO 490 L=1,200NSEP =0DO 500 I=3,NX+4DO 500 J=3,NY÷4IF (IMG(I,J).EQ.N) THENIMG(I,J) = N2DO 510 K=I-1,I+1,2IF (IMG(K,J).GT. 15) .AND. (IMG(K,J).LT. 31) THENIMG(K,J) N3NSEP = 1ENDIF510 CONTINUEDO 520 K=J-1,J÷1,2IF (IMG(I,K).GT. 15).AND.(IMG(I,K).LT.3 1) THENIMG(I,K) = N3NSEP = 1ENDIF520 CONTINUEENDIF500 CONTINUEIF (NSEP.EQ.0) GOTO 485490 CONTINUE485 CONTINUECC IN THE PREVIOUS LOOP, AFTER EACH PIXEL OF VALUE 10 OR 11C WAS CHECKED, IT WAS RE-SET TO 12 OR 14(10 - 12 AND 11-14)C SO THAT IT WOULD NOT BE CHECKED TWICE. NOW, THE PIXELSC WHICH WERE SET TO 12 AND 14 ARE RE-SET TO 10 AND 11.CDO 550 I=4,NX+3DO 550 J=4,NY+3IF (IMG(I,J).EQ.12) THENIMG(I,J) = 10ELSEIF (IMG(I,J).EQ.14) THEN141IMG(I,J) = 11ENDIFENDIF550 CONTINUECC PIXELS STILL SET TO 30 OR 15 ARE SET TO 10, THOSE SET TOC SEPRARESETTO11.CDO 240 I=3,NXi-4DO 240 J=3,NY+4IF (IMG(I,J).EQ.20) THENIMG(I,J) = 11ELSEif (IMG(I,J).EQ.30).OR.(IMG(I,J).EQ. 15) THENIMG(I,J) =10ENDIFENDIF240 CONTINUE360 CONTINUECwPJTE(9,*y’FINISHED SEPARATING PORES IN IMAGE “,NINCCC ALL PIXELS STILL SET TO 10 ARE RE-SET TO PORE VALUES,C AS ARE ANY PIXELS WITH VALUES GREATER THAN 30.C THOSE PIXELS SET TOll ARE SET TO SEPR (20)CDO 250 I=3,NX-t-4DO 250 J=3,NY÷4IF (IMG(I,J).EQ. 10).OR.(IMG(I,J).GT.30) THENIMG(I,J) = 255ELSEif (IMG(I,J).EQ.11) THENIMG(I,J) = 20ENDIFENDIF250 CONTINUECC THIS NEXT SECTION CHECKS TO ENSURE EACH PORE THROAT IS ATC LEAST TWO PIXELS IN WIDTH.CDO 270 M=1,2DO 260 I=4,NX+3DO 260 J=4,NY+3if (IMG(I,J).EQ.SEPR) THENNPOR=0IF (IMG(I+1,J).EQ.POR) NPOR=NPOR+1IF (JMG(I-l,J).EQ.POR) NPOR=NPOR+1IF (NPOR.EQ.2) THENIMG(I+1,J) = SEPRIMG(I- 1 ,J) = SEPRENDIFNPOR=0if (IMG(I,J÷1).EQ.POR) NPOR=NPOR+11 42IF (IMG(I,J-1).EQ.POR) NPOR=NPOR+1IF (NPOR.EQ.2) THENIMG(I,J÷1) = SEPRIMG(I,J-1) = SEPRENDIFENDIF260 CONTINUE270 CONTINUECC ALL PIXELS STILL SET TO 20 ARE ASSUMED TO BE WITHINC THE PORE THROATS AND ARE LEFT SET TO THIS VALUE.CCC THE SUBROUTINE PRO IS NOW CALLED TO LOCATE AND RETURNC THE LOCATIONS OF THOSE POINTS iN THE PORES WHICH AREC FARTHEST FROM THE EDGES OF THE PORES, JE THE CENTERSC OF THE PORES.C719 SEPR=20TMP = MDL(NINC)CALL PR(IMG,POR,TMP)WRITE(9,*)”AFTER PR”CC THE SUBROUTINE PRO IS NOW CALLED TO LOCATE AND RETURNC THE LOCATIONS OF THOSE POINTS IN THE PORE THROATS WHICH AREC FARTHEST FROM THE EDGES OF THE PORE THROATS, IE THE CENTERSC OF THE THROATS.CWRITE(9,*)”BEFORE PR”CALL PR(IMG,SEPR,TMP)WRITE(9,*)”AFTER PR”COPEN(6,FILE=OUTP(NINC))CC THE SUBROUTINE PRSZ IS CALLED TO DETERMINE THE SIZE OF THEC PORES AND PORE THROATS PRESENT IN THE IMAGE.CWRITE(9,*)’?BEFORE PRSZ”CALL PRSZ(IMG,PRDIST,THDIST,NPR,NTHR,MXRAD)WRITE(9,*)”AFrER PRSZ”CC THE SUBROUTINE SA IS CALLED TO DETERMINE THE POROSITY, POREC SURFACE TO VOLUME RATION, PORE THROAT SURFACE TO VOLUME RATIO,C AND THE SURFACE TO VOLUME RATIO OF THE TOTAL VOID SPACE.CCALL SA(IMG,POROS,PSV,TSV,THSV)wRrrE(9*)”AFTER SA”CC THE POROSITY, SURFACE TO VOLUME RATIOS, AND THE PORE AND POREC THROAT SIZES ARE WRHTEN INTO THE FILE PRSZ.”MDL”COUTP(NINC) = “PRSZ.”//MDL(NINC)OPEN(7,FILE=OUTP(NINC))143WR1TE(7, 109)109 FORMAT(”THE CURRENT IMAGE IS “,A20)WRITE(7,111) POROS111 FORMAT(”POROSITY OF MODEL “,F12.2)WRITE(7,1 12) PSV,THSV,TSV112 FORMAT(”SA/VOL RATIO FOR PORES ,F8 5/,”SA/VOL RATIO OF& PORE THROATS”,F8 5/,”SA/VOL RATIO OF TOTAL VOID SPACE tt,F8 5/)wRITE(7,*)”RAD PORE MN THROATWRITE(7,95)95 FORMAT(’ DIST DIST”,/i)DO 100 I=1,MXRADWRITE(7,96)I,PRDIST(I),THDIST(I)96 FORMAT(14,6X,15,4X,I5,4X,14)100 CONTINUECLOSE(7)CC THE IMAGE IS WRITTEN TO THE OUTPUT FILE.CDO 370 I=4,NX+3WRITE(6,9 1 )(IMG(I,J),J=4,NY÷3)370 CONTINUECC THE INPUT AND OUTPUT FILES ARE CLOSED.C13 CLOSE(S)CLOSE(6)CC THE NEXT STATEMENT IS THE END OF THE LOOP WHICHC IS REPEATED FOR EACH IMAGE.C380 CONTINUESTOPENDCC THE SUBROUTINES CALLED IN THE MAIN PROGRAM ARE INCLUDED.C PR.F CONTAINS PRO AND FINDS THE CENTERS OF PORES AND POREC THROATS, WRITES OUT THE CENTERS AND CHANGES THE PIXELC CORRESPONDING TO A CENTER TO THE VALUE 200.C SA.F CONTAINS SAO WHICH DETERMINES THE SURFACE TO VOLUMEC RATIO OF THE PORE SPACE AND THE POROSITY OF THE IMAGE.C PORSIZ.F CONTAINS PRSZQ AND DETERMINES THE MAX RADIUS OF AC CIRCLE WHICH IS CONTAINED ENTIRELY WITHIN THE PORE OR POREC THROAT AND CENTERED ON A PIXEL OF VALUE 200.C SPHTHR.F CONTAINS SPHTHR() IS CALLED BY PRSZO TO DO THEC COMPARISON BETWEEN THE IMAGE AND A REFERENCE CIRCLE GIVENC THE COORDINATES OF THE CENTER OF THE PORE OR PORE THROAT.C144INCLUDE PR.FINCLUDE SA.FINCLUDE PORSIZ.FINCLUDE SPHTHR.F145PRSZ.FCC MARCH, 1992 ALICE CHAPMANC ROCK PHYSICS LAB, UBCCC THIS SUBROUTINE ACCEPTS A CONFOCAL IMAGE FILE IN WHICHC THE PORES HAVE BEEN SEPARATED, AND THE PIXELS AT THEC CENTER OF PORES AND PORE THROATS HAVE BEEN RESET TOC THE VALUE 200. AN ARRAY CONTAINING A CIRCLE AND AC REFERENCE ARRAY ARE READ IN FROM THE FILE SPHS.CC THE CIRCLE IS PRODUCED BY SPHERE.F AND IS USED AS AC COMPARISON TO THE IMAGE TO DETERMINE THE SIZE OF THEC PORES OR PORE THROATS IN THE IMAGE GIVEN THE COORDINATESC OF THE POINT WITHIN THE PORES OR PORE THROATS WHICH ISC FARTHEST FROM THE EDGE.CCCC THE VARIABLES ARE DECLARED.CSUBROUTINE PRSZ(IMG,PRDIST,THDIST,NPR,NTHR,MXRAD)CINTEGER*2 IMG(530,780), SPH(200,200)INTEGER SOL,SEPR,RAD,RAD2,CTR,MXRADINTEGER PRDIST( 101),THDIST(101)INTEGER NTHR,CHKSPH( 101) ,NPRREAL RFIT,MXFITCHARACTER*20 MDL,INPCC***************************CC VARIABLE DESCRIPTION:CC IMG(530,780): THE CONFOCAL IMAGE IS PLACED IN THIS ARRAYC WITH THREE ROWS OF SOLiDS SURROUNDING ITC SPH(200,200): AN ARRAY CONTAINING A SET OF CONCENTRICC CIRCLES WITH THE OUTSIDE ROW OF PIXELS IN EACHC CIRCLE BEING NUMBERED (RAD+0.5) WERE RAD WAS THEC RADIUS OF THE CIRCLE IN QUESTION. IE: THEC CENTERMOST PIXEL IS NUMBERED 1 SINCE RAD = 0.5.C ALSO, THE CENTER IS AT SPH(100.5,100.5)C CHKSPH(I): A COMPARISON ARRAY WHICH RECORDS HOW MANYC PIXELS RE WITHIN A CIRCLE OF RADIUS (1-0.5)C SOL: THE NUMBER WHICH PIXELS WITHIN THE SOLID PORTIONC OFTHEIMAGEWILLBESETTO(0)C POR: THE NUMBER WHICH PIXELS WITHIN THE PORES OF THEC IMAGE WILL BE SET TO (255)C SEPR: THE VALUE WHICH PIXELS IN THE PORE THROATS WILLC BE SET TO TO SEPARATE THE PORES (20)C CrR: THE NUMBER TO WHICH THE PIXEL AT THE CENTRE OF AC PORE OR PORE THROAT IS SETC NPR: THE NUMBER OF PORES WHOSE SIZE WAS MEASURED146C NTHR: TH NUMBER OF PORE THROATS WHOSE SIZE WAS MEASUREDC RAD: THE SIZE OF THE CIRCLE WHICH FITS ENTIRELY WITHINC THE PORE OR THROAT BEING MEASURED. (SIZE: ACTUALC RADIUS +0.5)C PRDISTjTHDIST: AFTER THE SIZE OF THE PORE OR PORE THROATC IS DETERMINED THE SIZE IS RECORDED IN PRDIST (PORES)C AND THDIST (PORE THROATS). PRDIST(RAD) OR THDIST(RAD)C IS INCREMENTED BY ONE AFTER EACH MEASUREMENT.C MXRAD: THE INTEGER REPRESENTING THE SIZE OF THE LARGEST POREC SIZE MEASURED. MXRAD = RAD1US-i-0.5C MDL: THE NAME OF THE MODEL, WHICH IS USED TO NAME SOMEC OUTPUT FILESCC***************************CCC THE WIDTH (NX) AND LENGTH (NY) OF THE CONFOCAL IMAGE ARRAYC AND THE WORKING ARRAY (NTX) AND (NTY) ARE SET.CNX = 512NY = 768NTX = NX + 6NTY = NY +6CC THE ARRAYS IN THE PROGRAM ARE INiTIALIZED TO ZERO.CDO 10 1=1,101PRDIST(I) =0THDIST(I) 0CHKSPH(I) = 010 CONTINUECC THE SOLID, POR, SEPARATION, AND PORETHROAT CENTER VALUES FORC THE IMAGE FILE ARE INITIALIZED.CSOL=0SEPR =20POR = 255CTR = 200CC THE COMPARISON SPHERE AND ARRAY (CHKSPH) ARE NOW READ IN FROMC THEFILE”SPHS’CINP = ‘SPHS’OPEN(7,FILE=INP)DO 201=1,200READ(7,91 )(SPH(I,J),J=1 ,200)20 CONTINUE91 FORMAT(10014)CDO 301=1,101READ(7,94)J,CHKSPH(I)IF (I.NE.J) THEN147WPJTE(9*) “THERE IS A PROBLEM WITH READING IN THE”WRITE(9,*)”CHKSPH ARRAY FROM THE FILE: SPHS.”WRITE(9,*)”DO YOU WISH TO QUIT? IF SO, ENTER A 1”READ(9,*)KIF (K.EQ.1) THENCLOSE(7)GOTO 11ENDIFENDIF30 CONTINUE94 FORMAT(218)CCLOSE(7)CC THE COUNTING VARIABLES USED TO KEEP TRACK OF THE MAXIMUMC RADIUS OF AND THE NUMBER OF PORE AND PORE THROATS MEASUREDC ARE INITIALIZED TO ZERO.CMXRAD =0NPR=0NTHR =0CC THE SIZE OF THE PORES AND PORE THROATS IN WHICH THE POINTSC FARTHEST FROM THE EDGES OF THE PORES HAVE BEEN SET TO CTRC ARE NOW MEASURED.CDO 40 I=4,NX+3DO 40 12=4,NY+3IF (IMG(I,12).EQ.200) THENJ=IK =12CC SPHTHR IS CALLED TO DETERMINE THE SIZE OF THE LARGEST CIRCLEC CENTERED AT IMG(J,K) WHICH FITS ENTIRELY WITHIN THE PORE/THROAT.C THE SIZE WHICH IS RETURNED IS CHECKED AND, IF IT IS LARGER THANC ANY PREVIOUSLY RETURNED, MXRAD IS RE-SET TO EQUAL RAD.CCALL SPHTHR(J,K,IMG,SPH,NREF,RFIT,CHKSPH)IMG(I,12) = 210RAD = INT(RFIT)IF (RAD.GT.MXRAD) MXRAD = RADIF (RAD.GT.0) THENCC IF THE REFERENCE NUMBER RETURNED BY SPHTHR IS EQUAL TO POR,C THEN THE SIZE OF A PORE WAS MEASURED AND NPR AND PRDIST(RAD) AREC INCREMENTED. IF NREF WAS EQUAL TO SEPR, A PORE THROAT WASC MEASURED AND NTHR AND THDIST(RAD) ARE INCREMENTED. IF NREF WASC 0, THOUGH, THIS INDICATES THE PORE OR PORE THROAT WAS IN CONTACTC WITH THE EDGE OF THE IMAGE AND THEREFORE ITS MEASUREMENT WILLC NOT BE RECORDED.CIF (NREF.EQ.POR) THENNPR=NPR+ 1PRDIST(RAD) = PRDIST(RAD) + 1148ELSEIF (NREF.EQ.SEPR) THENNTHR = NTHR + 1THDIST(RAD) = THDIST(RAD) + 1ENDIFENDIFEND1FENDIF40 CONTINUECC THE IMAGE IS NOW CHECKED FOR PORE/PORE THROATS THAT WEREC MEASURED BUT WHOSE CENTERS WERE LATER INCLUDED IN ANOTHERC PORE AND THEREFORE THEMEASUREMENT MUST BE REMOVED FROM THEC DISTRIBUTION.CDO 50 I=4,NX÷3DO 50 12=4,NY+3IF (IMG(I,I2).EQ.220) THENJ=IK=12CC SPHTHR IS CALLED TO DETERMINE THE SIZE OF THE LARGEST CIRCLEC CENTERED AT IMG(J,K) WHICH FITS ENTIRELY WITHIN THE PORE/THROAT.CCALL SPHTHR(J,K,IMG,SPH,NREF,RFIT,CHKSPH)IMG(I,I2) = NREFIF (NREF.EQ.0) THENIMG(I,I2) = IMG(I+1,12)ENDIFRAD = INT(RFIT)IF (RAD.GT.0) THENCC IF THE REFERENCE NUMBER (NREF) RETURNED BY SPHTHR IS EQUAL TOC POR, THEN THE SIZE OF A PORE WAS MEASURED AN]) NPR ANDC PRDIST(RAD) ARE DECREASED BY 1. IF NREF WAS EQUAL TO SEPR, A POREC THROAT WAS MEASURED AND NTHR AN]) THDIST(RAD) ARE DECREASED BYC 1. IF NREF WAS 0, THOUGH, THIS INDICATES THE PORE OR PORE THROATC WAS IN CONTACT WITH THE EDGE OF THE IMAGE AND THEREFORE ITSC MEASUREMENT WILL NOT BE RECORDED.CIF (NREF.EQ.POR) THENNPR=NPR- 1PRDIST(RAD) = PRDIST(RAD) - 1ELSEIF (NREF.EQSEPR) THENNTHR = NTHR -1THDIST(RAD) = THDIST(RAD) - 1ENDIFENDIFENDIFELSEIF (IMG(I,12).EQ.210) THENIMG(I,12) = 200ELSE149IF (IMG(I,12).EQ.205) THENIF (IMG(I-1,12).GT.O) THENNREF = IMG(I- 1,12)ELSEIF (IMG(I+1,12).GT.O) THENNREF = IMG(I+1,12)ENDIFENDIFIMG(I,12) NREFENDIFENDIFENDIF50 CONTINUE11 RETURNEND150PR.FCC MARCH, 1992 ALICE CHAPMANC ROCK PHYSICS LAB, UBCCC THIS SUBROUTINE TAKES CONFOCAL IMAGES IN WHICHC THE PORES HAVE BEEN SEPARATED AND FINDS THE PIXELSC WiTHIN THE PORES OR PORE THROATS WHICH ARE FARTHESTC FROM THE EDGES. IF THE VALUE OF NREF IS 255, THEC VALUE CORRESPONDING TO PORES N THE IMAGE, THE CENTERC OFTHEPORESAREFOUNDAND,IFITIS2O,THECENTEROFC THE PORE THROATS ARE FOUND.CC THE POINT IN THE CENTERMOST PORTION OF THE PORES ANDC PORE THROATS IS FOUND USING THE WATERSHED ALGORITHM.C N THIS ALGORITHM, THE PIXELS ON THE OUTER EDGES OFC THE PORE/THROAT ARE SET TO 50, THEN INNER ONES IN CONTACTC WITH THE 50’S ARE SET TO 51...ETC., UNTIL THE PORES/THROATSC ARE COMPLETELY RE-NUMBERED. THEN THOSE PIXELS GREATERC THAN 50 WHICH DO NOT HAVE A LARGER NUMBER WITHIN A BOXC WHICH EXTENDS 5 PIXELS ON EACH SIDE OF THEM ARE ASSUMEDC TO BE NEAR THE CENTERS. THE PIXELS UP TO 20 PIXELS ONC EITHER SIDE OF SUCH A PIXEL ARE CHECKED TO DETERMINEC IF THEY ARE PART OF A CLUSTER OF SIMILARLY VALUED PIXELSC INCLUDiNG THE ORIGINAL PIXEL. IF SO, THE COORDINATES OFC THE PIXELS IN THE CLUSTER ARE AVERAGED TO DETERMINE THEC ACTUAL CENTER OF THE PORE/THROAT.CCCCSUBROUTINE PR(IMG,NREF,MDL)CC THE VARIABLES ARE DECLARED.CINTEGER*2 IMG(530,780)INTEGER POR,SOL,SEPR,CTR( 10000,2),BDRYINTEGER NSPH,MAX,NREFCHARACTER*20 MDL,OUTPCC***************************CC VARIABLE DESCRIPTION:CC IMG(530,780): THE CONFOCAL IMAGE IS PLACED IN THIS ARRAYC WITH THREE ROWS OF SOLIDS SURROUNDING ITC POR: THE NUMBER WElCH PIXELS WITHIN THE PORES OF THEC IMAGE ARE SET TO (255)C SOL: THE NUMBER WHICH PIXELS WITHIN THE SOLID PORTION OFC THE IMAGE ARE SET TO (0)C SEPR: THE NUMBER WHICH PIXELS IN THE SMALLEST PART OF THEC PORE THROATS ARE SET TO TO SEPARATE THE PORES (20)C NREF:THE NUMBER WHICH INDICATES WHETHER PORE OR PORE THROAT151C CENTERS WILL BE DETERMINED. IF NREF = POR, THEN POREC CENTERS WILL BE FOUND, AND IF NREF = SEPR, PORE THROATC CENTERS WILL BE FOUND.C CTR(K,I): A WORKING ARRAY WHICH HOLDS THE COORDINATES OF THEC PIXELS WiTHIN THE PORES OR PORE THROATS WHICH AREC FARTHEST FROM THE EDGES OF THE PORES.C NSPH:THE TOTAL NUMBER OF PORES OR THROATS WHOSE CENTERS AREC RECORDED IN CTR.C BDRY:THE INfl1AL NUMBER WHICH PORE PIXELS BORDERING THE SOLIDC ARE SET TO IN THE WATERSHED ALGORITH1vI. PORE PIXELS WHICHC BORDER THOSE WHICH ARE SET TO BDRY ARE THEN SET TO BDRY-i-1C ETC, ETC.C MAX: THE LARGEST NUMBER USED IN THE WATERSHED ALGORITHM. (IE:C THE LARGEST NUMBER WHICH A PORE PIXEL IS RE-SET TO.C MDL: THE NAME OF THE MODELS, USED TO NAME OUTPUT FILESCC***************************CC THE WIDTH (NX) AND LENGTH (NY) OF THE CONFOCAL IMAGE ARRAYC AND THE WORKING ARRAY (NTX) AND (NTY) ARE SET.CNX = 512NY = 768NTX = NX + 6NTY = NY +6CC THE ARRAYS IN THE PROGRAM ARE INITIALIZED TO ZERO.CDO 201=1,2DO 20 J=1,10000CTR(J,I) = 020 CONTINUECC THE PORE/SOLID, SEPARATION, AND BOUNDARY VALUES FOR THE IMAGEC FILE ARE INITIALIZED.CPOR=255SOL=0SEPR =20BDRY =50CC ThE “WATERSHED ALGORiTHM IS NOW APPLIED. THIS SETS THE VALUESC OF THE PIXELS WITHIN THE PORES TO PROGRESSWELY LOWER OR HIGHERC VALUES DEPENDING ON THEIR PROXIMITY TO THE SOLID.CWRITE(9,*)”STARTING WTSHD PORTION”DO 40 I=BDRY,255NTST=0ITMP = IDO 50 J = 4,NX+3DO 50 K= 4,NY+3IF (I.EQ.BDRY) THENIF (IMG(J,K).EQ.POR).OR.(IMG(J,K).LT.25) THENIF (IMG(J,K).NE.NREF) THEN1 52NTST= 1IF (IMG(J+1,K).EQ.NREF) IMG(J+1,K) = ITMPiF (IMG(J-1,K).EQ.NREF) IMG(J-1,K) = ITMPIF (IMG(J,K+1).EQ.NREF) IMG(J,K+1) = ITMPIF (IMG(J,K-1).EQ.NREF) IMG(J,K-1) = ITMPENDIFENDIFELSEIF (IMG(J,K).EQ.(ITMP-1)) THENNTST= 1IF (IMG(J+1,K).EQ.NREF) IMG(J+1,K) = ITMPIF (IMG(J-1,K).EQ.NREF) IMG(J-1,K) = ITMPIF (IMG(J,K+1).EQ.NREF) IMG(J,K+1) = ITMPIF (IMG(J,K-1).EQ.NREF) IMG(J,K-1) = ITMPENDIFENDIF50 CONTINUEIF (NTST.EQ.0) THENMAX = I-iGOTO 111ENDIF40 CONTINUECC THE POSSIBLE CENTERS OF PORES ARE NOW SELECTED, THEIRC COORDINATES STORED IN THE ARRAY CRC111 NSPH=0WPJTE(9,*)”PAST WTSHD PORTION”DO 751= MAX,BDRY,-1DO 70 J=4,NXi-3DO 70 K=4,NY+3CC EACH PIXEL IN THE IMAGE IS CHECKED UNTIL ONE WHICH HASC BEEN ALTERED BY THE WATERSHED ALGORITHM IS FOUND. THEC PARAMETER NPOR IS THEN SET TO THE VALUE OF THIS PIXELC AND THE TEST VARIABLE, 1TST IS SET TO 0.CIF (IMG(J,K).EQ.I) THENNPOR = IMG(J,K)ITST=0CC THIS LOOP CHECKS A BOX IN THE IMAGE WHICH EXTENDS 5 PIXELSC TO EACH SIDE OF THE PIXEL SELECTED. IF THE VARIABLES J2,C K2 GO OUT OF RANGE, THE PIXEL IS NOT CHECKED, OTHERWISE iTC IS CHECKED TO DETERMINE IS IT IS LARGER THAN THE CHOSENC PIXEL’S VALUE OF NPOR. IF SO, ITST IS INCREMENTED BY 1.CDO 80 J2=J-5,J-t-5DO 80 K2=K-5,K+5IF ((J2.LT.l).OR.(J2.GT.NTX)) GOTO 80IF ((K2.LT.1).OR.(K2.GT.NTY)) GOTO 80IF (IMG(J2,K2).GT.NPOR) THEN1 53ITST = ITST+1ENDIF80 CONTINUECC IF THERE WERE NO PIXELS LARGER THAN THE ONE SELECTED INC THEAREAAROUNDIT,THENITISASSUMEDTOBEATORNEARC A PORE/PORE THROAT CENTER AND ThE NUMBER OF CENTERS FOUNDC , NSPH, IS INCREMENTED BY ONE. THE POSSIBLE CENTER, ORC ORIGIN OF THE PORE/THROAT IS THEN SET TO THE COORDINATESC OF THE CHOSEN PIXEL. (ORU = J, ORIK = K) NORI IS THEC NUMBER OF PIXEL POSITIONS AVERAGED TO DETERMINE THE ACTUALC CENTER OF THE PORE THROAT AND IS INiTIALIZED TOP 1.CIF (ITST.EQ.0) THENNSPH = NSPH + 1ORIJ=JORIK=KNORI= 1CC THE AREA AROUND THE PIXEL IN A BOX EXTENDING UP TO 20 PIXELSC TO EITHER SIDE IS NOW CHECKED AND ANY PIXELS OF THE SAMEC VALUE ARE AVERAGED TO DETERMINE THE ACTUAL CENTER OF THEC PORE/THROAT. ONCE THE PIXEL HAS BEEN ALTERED, iT IS RE-SET TOC 240 SO THAT ThE SAME PORE/THROAT CENTER WILL NOT BE “FOUND” FROMC MORE THAN ONE ORIGINATING PIXEL.C FOR EACH LOOP IN WHICH L IS INCREMENTED, THE PORTION OF THEC BOX TO BE CHECKED IS ONE PIXEL FARTHER FROM THE CENTER THANC BEFORE AND WITHIN EACH LOOP, ITST IS ORIGINALLY SET TOO, ANDC IF A PIXEL OF VALUE NPOR IS FOUND, RE-SET TO 1. IF NO PIXELSC OF VALUE NPOR ARE FOUND, THE LOOP IS TERMINATED.CDO 90 L=1,20N2=L*21TST=0DO 100 J2=J-L,J+L,N2DO 100 K2=K-L,K+LCC THE NEXT TWO STATEMENTS CHECK IF J2 OR K2 IS OUT OF RANGEC OF THE ARRAY BOUNDARIES AND, IF SO, SKIPS CHECKING THEC VALUES OF THE PIXELS AT THOSE POSITIONS.CIF ((J2.LT.2).OR.(J2.GT.NTX-1)) GOTO 100IF ((K2.LT.2).OR.(K2.GT.NTY-1)) GOTO 100IF (IMG(J2,K2).EQ.NPOR) THENDO 105 J3=J2-1,J2+1DO 105 K3 = K2-1,K2+1IF (IMG(J3,K3).EQ.240) THENORIJ = (ORIJ*REAL(NORI) + REAL(J2))/(1.O -i-REAL(NORI))ORIX = (ORIK*REAL(NORI) + REAL(K2))/(1.O +REAL(NORI))NORI=NORI÷ 1ITST= 1GOTO 100ENDIF105 CONTINUE1 54ENDIF100 CONTINUEDO 110 K2=K-L,K+L,N2DO 110 J2=J-L+1 ,J+L- 1CC THE NEXT TWO STATEMENTS CHECK IF J2 OR K2 IS OUT OF RANGEC OF THE ARRAY BOUNDARIES AND, IF SO, SKIPS CHECKING THEC VALUES OF THE PIXELS AT THOSE POSITIONS.CIF ((J2.LT.2).OR.(J2.GT.NTX-1)) GOTO 110IF ((K2.LT.2).OR.(K2.GT.NTY-1)) GOTO 110IF (IMG(J2,K2).EQ.NPOR) THENDO 115 J3=J2-1,J2+1DO 115 K3=K2-1,K2-f-1IF (IMG(J3,K3).EQ.240) THENORIJ = (ORIJ*REAL(NORI) + REAL(J2))/(1.0 -i-REAL(NORI))ORIK = (ORIK*REAL(NORI) + REAL(K2))/(1.0 +REAL(NORI))NORI=NORI+ 1ITST = 1IMG(J2,K2) = 240GOTO 110ENDIF115 CONTINUEENDIF110 CONTINUECC IF ITST IS STILL EQUAL TOO, THEN NO PIXELS OF VALE NPOR HAVEC BEEN LOCATED IN THIS LOOP AND ALL PIXELS SURROUNDING THE POREC CENTER ARE ASSUMED TO HAVE BEEN FOUND SO THE LOOP ISC TERMINATED.CIF (ITST.EQ.0) GOTO 12190 CONTINUECC THE J AND K COORDINATES OF THE CENTER OF THE PORE/THROAT AREC ROUNDED UP IF THE DECIMAL PORTION IS GREATER THAN 0.5, DOWNC OTHERWISE AND THEN CHANGED TO INTEGERS AND PLACED IN THE CTRC ARRAY.C121 IF ((ORIJ + O.5).GT.(INT(ORIJ+O.5))) THENORU = ORIJ + 0.5ENDIFIF ((ORIK + O.5).GT.(INT(ORIK+0.5))) THENORIK=ORIK+0.5ENDIFCrR(NSPH,1) = INT(ORIJ)CTR(NSPH,2) = INT(ORIK)IMG(INT(ORIJ),INT(ORIK)) = 200ENDIFENDIF70 CONTINUE75 CONTINUEWRITE(9,*)“FINISHED WITH CENTERS”C1 55C THE ARRAY CTR IS NOW CHECKED TO ENSURE THAT REPEATED CTRC COORDINATES ARE REMOVED AND THAT ANY COORDiNATES WHICHC CORRESPOND TO PORE THROATS ARE ALSO REMOVED.CDO 160 L= 1,NSPHIF (CTR(L,1).GT.0) THENIF (CTR(L,2).GT.0) THENM = CTR(L,1)N = CTR(L,2)DO 170 J=L÷1,NSPHIF (CTR(J,1).EQ.M) THENIF (CTR(J,2).EQ.N) THENCTR(J,1) =0CTR(J,2) =0ENDIFENDIF170 CONTiNUEENDIFENDIF160 CONTINUECC THE ARRAY CTR IS NOW CHECKED TO ENSURE THAT CTR COORDINATESC WHICH ARE NOT VIABLE (IE: LT 0 OR EQUAL TO ANOTHER CUR) AREC REMOVED.CNCHK=0DO 180 I=1,NSPHIF (CTR(I,1).GT.0) THENIF (CTR(I,2).GT.0) THENNCHK = NCHK + 1IF (NCHK.LT.I) THENCTR(NCHK,1) = CTR(I,1)CTR(NCHK,2) = CTR(I,2)CTR(I,1) =0CTR(I,2) = 0ENDIFENDIFENDIF180 CONTINUENSPH = NCHKCC THE PIXELS WHICH HAVE NOT BEEN RE-SET TO THEIR ORIGINALC VALUE OF NREF ARE NOW RE-SET.CDO 300 I=4,NX+3DO 300 J=4,NY+3IF (IMG(I,J).GT.SEPR).AND.(IMG(I,J).NE.200) THENIF (IMG(I,J).NE.POR) THENIMG(I,J) = NREFENDIFENDIF300 CONTINUE1 56CC THE PIXELS WHICH CORRESPOND TO THE CENThR COORDINATESC DETERMINED WITH THIS PROGRAM ARE NOW SET TO THE VALUEC 200CDO 310 I=1,NSPHJ = CTR(I,1)K = CTR(I,2)IMG(J,K) = 200310 CONTINUERETURNEND1 57SA.FCC MARCH, 1992 ALICE CHAPMANC ROCK PHYSICS LAB, UBCCC THIS SUBROUTINE ACCEPTS A CONFOCAL IMAGE FILE iN WHICH THEC PORES HAVE BEEN SEPARATED. THE SURFACE AREA, VOLUME,C AND THE SURFACE AREA TO VOLUME RATIO OF THE ENTIRE IMAGEC IS DETERMINED AS WELL AS THE POROSITY.CCC***************************CSUBROUTINE SA(IMG,POROS,PSV,TSV,THSV)CC THE VARIABLES ARE DECLARED.CINTEGER*2 JMG(530,780)INTEGER SOL,SEPR,TVOL,TSA,PVOL,THVOLINTEGER PSA,THSAREAL POROS ,PSV,TSV,THSVCHARACTER*20 MDLCC***************************CC VARIABLE DESCRIPTION:CC IMG(530,780): THE CONFOCAL IMAGE IS PLACED IN THIS ARRAYC WITH THREE ROWS OF SOLIDS SURROUNDING ITC SOL: THE NUMBER WHICH PIXELS WITHIN THE SOLD) PORTION OFC THEIMAGEWILLBESETTO(0)C SEPR: THE NUMBER WHICH PIXELS IN THE SMALLEST PART OF THEC PORE THROATS WILL BE SET TO TO SEPARATE THE PORES (20)C POR: THE NUMBER WHICH PIXELS WITHIN THE PORES OF THE IMAGEC WILL BE SET TO (255)C CrR: THE NUMBER TO WHICH THE PIXEL AT THE CENTRE OF A POREC OR PORE THROAT IS SETC PVOL: THE NUMBER (VOLUME) OF PORE PIXELS IN THE IMAGEC THVOL:THE NUMBER (VOLUME) OF PORE THROAT PIXELS IN THE IMAGEC TVOL: THE TOTAL NUMBER (VOLUME) OF PIXELS IN THE IMAGEC PSA: THE SURFACE AREA OF PORES IN THE IMAGEC THSA: THE SURFACE AREA OF PORE THROATS IN THE IMAGEC TSA: THE TOTAL SURFACE AREA OF PORES AND THROATS IN THEC POROS: THE POROSITY OF THE IMAGE: IE THE RATIO OF THE NUMBERC OF iMAGE PORE AND PORE THROAT PIXELS IN THE IMAGE TOC THE TOTAL NUMBER OF PIXELS IN THE IMAGE.C PSV: THE SURFACE TO VOLUME RATIO (PSA/PVOL) OF THE PORESC THSV: THE SURFACE TO VOLUME RATIO (THSAITHVOL) OF THE POREC THROATSC TSV: THE OVERALL SURFACE TO VOLUME RATIO OF THE VOID SPACEC IN THE IMAGE: (PSA + THSA)/(PVOL + THVOL)C MDL: THE NAME OF THE MODEL, WHICH IS USED TO NAME SOMEC OUTPUT FILESC1 58c***************************CC THE WIDTH (NX) AND LENGTH (NY) OF THE CONFOCAL IMAGE ARRAYC AND THE WORKiNG ARRAY (NTX) AND (NTY) ARE SET.CNX = 512NY = 768NTX = NX + 6NTY = NY +6CCC THE SOLID, SEPARATION, PORE AND CENTER VALUES FOR THE IMAGE FILEC ARE INITIALIZED.CSOL=0SEPR =20POR=255CTR=200CC THE COUNTERS FOR THE VOLUMES AND SURFACE AREAS OF PORES AND POREC ThROATS IN THE MODEL ARE IN111ALIZED TO ZERO.CTVOL=0PVOL=0THVOL=0PSA =0THSA=0TSA=0CC THE FOLLOWING LOOP CHECKS EACH PIXEL IN THE IMAGE AND INCREMENTSC THE VOLUMES AND SURFACE AREA COUNTERS APPROPRIATELY.CDO 10 I=4,NX+3DO 10 J=4,NY+3CC THE TOTAL VOLUME COUNTER (TVOL) IS INCREMENTED FOR EACH PIXEL.CTVOL = TVOL + 1ITST=0CC THIS SERIES OF IF/ThEN STATEMENTS SETS ITST TO 2 IF THE PIXEL ISC IN A PORE AND TO 3 IF iT IS IN A PORE THROAT.CIF (IMG(I,J).EQ.CTR) THENDO 20 12=1-1,1+1DO 20 J2=J- 1 ,J+ 1IF (IMG(I,J).EQ.POR) THEN1TST=2ELSEIF (IMG(I,J).EQ.SEPR) 1TST =3ENDIF20 CONTINUEELSEIF (IMG(I,J).EQ.POR) THEN1 59ITST =2ELSEIF (IMG(I,J).EQ.SEPR) 1TST = 3ENDIFENDIFCC IFTHEPIXELWASINAPORETHROAT,ITSTIS3ANDTHVOLISC INCREMENTED BY ONE. THE FOUR ADJOINING PIXELS ARE THEN CHECKEDC AND, IF THE PIXEL CHECKED IS A SOLID, THSA IS INCREMENTED BY 1.C SIMILARLY, IT THE PIXEL WAS IN A PORE, ITST IS 2, AND PVOL ISC INCREMENTED BY ONE. THE FOUR ADJOINING PIXELS ARE THEN CHECKEDC AND, IF THE PIXEL CHECKED IS A SOLID, PSA IS INCREMENTED BY 1.CIF (ITST.EQ.3) THENTHVOL = THVOL + 1IF (IMG(I÷1,J).EQ.SOL) THSA = THSA + 1IF (IMG(I-1,J).EQ.SOL) THSA = THSA + 1IF (IMG(I,J+1).EQ.SOL) THSA = THSA + 1IF (IMG(I,J-1).EQ.SOL) THSA = THSA + 1ELSEIF (ITST.EQ.2) THENPVOL = PVOL + 1IF (IMG(I+1,J).EQ.SOL) PSA = PSA + 1IF (IMG(I-1,J).EQ.SOL) PSA = PSA + 1IF (IMG(I,J+1).EQ.SOL) PSA = PSA + 1IF (IMG(I,J-1).EQ.SOL) PSA = PSA + 1ENDIFENDIF10 CONTINUECC THE TOTAL SURFACE AREA (TSA) OF THE VOID SPACE IS CALCULATED.CTSA = PSA + THSACC THE POROSITY AND THE SURFACE AREA TO VOLUME RATIO OF THE PORESC AND PORE THROATS ARE CALCULATED AND PRINTED OUT.CTSV = REAL(TSA)/REAL(PVOL + THVOL)PSV REAL(PSA)/REAL(PVOL)THSV = REAL(THSA)/REAL(THVOL)POROS = REAL(PVOL + ThVOL)/REALfVOL)RETURNEND1 601.91.333333333333333:cj33333333333333:N0IJpJ3SGIHVfl1VA1HivanflN’DIN‘TfN’HN’(loJ)HJSN”TcI’Nu‘XVNJI9aLNII1N‘(JoT)HaSNHD‘1SL‘GVIITIDLLNI(ooVoo)Has‘(o8L’oc)ov’Izrnou’..ii(HasiH3’1IiI’iaIN’Has’oNI’IIN’IIN)llUHdSItll1flO1HflS******************************V13iDHfl‘fiVISDISAHJ)J3(j{:yJAJ1}{333J’prS1V011iL310d103ZIS31113MIELL3(I01A1TlVI’ 41ThdC135flSISITUIfN’TINSLLVNJcflI0ODTc1UN3DLDNISflS3I0d3HLMIG35013N3KITWrIcIWO3SIHDIHML313IIDIJiLI03ZIS3H1S3M1WHL3GN0IS’I3AIV’1fl3LLIVdSITUN0LLSflONIS1IXIJ3HL3C[fflDNI(flflOMH3fflA313II3IIU10co÷sniciviKLoiGNdSIOD5I3mNnN3111•S3oI3ZioaNnoI0N3vaVNOSI3HNflNMI(ULLV3IGMI3IVHJSNIS3’T3II3HL1V0IIU3I0cT/3I0J311101HJSXIaLVN3111NIRLIMU31’IdcIflSS313II3‘EIHLS3IVJNO3‘1V0IHL3I0dI03110dVtIO‘IaLM333111I0SLVMCII0ODAQNVX3H1N3A19N31IA‘TllALflOIHflS/NVIO0JcISIRL(co-i)SflIUVN10fl3ll3VN11UJAkaIV1V0IRLflI0dI0310J 311110S’13X1JANVNAkOHSCfI033NHDIHAkAVfflVNOSIIIVdWO3y:(I)HdSN9J411‘AVINVEI9VVJI31{LNI1V0llU3I0J10aI0d3H1dOITLLN3C)1RL10S1VNK0033HJ,:TIN’UKSGMflOHAVRIV3H1NllUI/*ATILLSSIH3IHMITLLN3331-ILIA10H3DMVISKI1S3H1IVd311105WSIllSflHI‘HdSNI313TID3DN3Ia1aI1S30aVl3H1dOSflICIVI3RLSII+XVN:yiHCTSAV’RIV3111NIS313ll333N3lad3I3111dOIaLM333111dOSLVNIl003:TTT)TJJ1V0ll1L3l0dVSVM1J’O=daINllGNV3I0dVSVM.IIN3H1‘cc=d3INdlU3IflSV3NM33HSVH1V0lHL3I0dI03I0JVI3RL3HMSLV3KIMI3fflAfl43fTJ(co+snicivi1VfLL3V:3ZIS)UlflSV3N9N13H1V0IHII03l0d3111MIHIIMA13ULN3SlIdHDIHAkWT3llD3111dO3ZIS3111(co-i)SflIUVINIV1I33VdO3’1DllDVNIIUIM3H(flflOHSS13XIJANVNA&OHS([I0D3IHDIIIMAVIVN0SNVcWI0DV:(j)j(cool’cool)HaSIVSIIaLN3D3111co=(IVI3DNISI0113551‘13X1J1SowIaLM3D3HV31NOIIS3flONI313II3RLdOSfllCIV13H1SIC[VI3l3HhX(co+uva)(J3i3HWflNONI3H3’I3ll3HDV3NI513X1JdOM0l3GISIflO3111HJJh.S3’IDII3DIULM3DNODdO13SV9NJNIVIMODAVIIVNVIIONIc[MflOIflS501105dOSM0I331RLHJJMAVR1VSITUNI(133VicISI3OVIAII1VDOdNOD3111:(O8L’oEc)9jAJ333333333333333C***************************CC THE CENTER COORDINATES OF THE REFERENCE CIRCLE (NK1,NL1), THEC MAXIMUM DIMENSIONS OF THE IMAGE ARRAY (MXX,MXY), AND THEDISTANCEC FROM THE CENTER OF THE REFERENCE CIRCLE TO THE EDGE OF THEARRAYC SPH(MAX) ARE INiTIALIZED.CNK1 = 100NL1 = 100RFIT = 1MXFIT= 1MAX =99MXX = 530MXY=780CC THE ARRAY NSPHQ IS INITIALIZEDCDO 10 1=1,101NSPH(I) = 010 CONTINUECC THE AREA IN A BOX EXTENDING UP TO 20 PIXELS TO EiTHER SIDE OF THEC CENTER IF NOW CHECKED AND PIXELS WHICH ARE PORE OR PORE THROATC PIXELS RESULT IN THE APPROPRIATE ARRAY (NSPH) BEING INCREMENTEDC BY ONE. THE ARRAY POSITION WHICH IS INCREMENTED IS THAT WHICHC CORRESPONDS TO THE RADIIJS+1.5 OF THE LARGES CIRCLE CONTAININGC THE PIXEL WHICH WAS CHECKED.C FOR EACH LOOP IN WHICH L IS INCREMENTED, THE PORTION OF THE BOXC TO BE CHECKED IS ONE PIXEL FARTHER FROM THE CENTER THAN IN THEC PREVIOUS LOOP. WITHIN EACH LOOP, NSPH(L+2), THE MAXIMUM SIZE OFC CIRCLE WHICH WOULD FIT WITHIN THE AREA CHECKED IN THAT LOOP, ISC COMPARED TO THE REFERENCE ARRAY CHKSPH(L+2) AND, IF THEY AREC NOT EQUAL, THE LOOP IS TERMINATED.CDO2OL= 1,MAXN2 = L*2NSZ=L+2DO 30 NI = NI1-L,M1-f-L,N2IF ((NI.LT.1).OR.(NLGT.MXX)) GOTO 30DK = NI-MiNK=NK1 +DKDO 30 NJ = NJ1-L,NJ1+LIF ((NJ.LT.1).OR.(NJ.GT.MXY)) GOTO 30DL=NJ-NJ1NL=NL1 +DLNIN = SPH(NK,NL) + 1IF (IMG(NI,NJ).EQ.200) THENIMG(M,NJ) = 205ELSEIF (IMG(NI,NJ).EQ.210) THENIMG(NI,NJ) = 220ENDIF1 62ENDIFIF (IMG(NI,NJ).GT.0) THENNSPH(NIN) = NSPH(NIN) + 1ENDIF30 CONTINUEDO 40 NI = NI1-L+1,NI1+L-1IF ((NI.LT.1).OR.(NI.GT.MXX)) GOTO 40DK = NI-NilNK=NK1 +DKDO 40 NJ = NJ1-L,NJ1+L,N2IF ((NJ.LT.1).OR.(NJ.GT.MXY)) GOTO 40DL=NJ-NJ1NL=NL1 +DLIF (IMG(NI,NJ).EQ.200) THENIMG(NI,NJ) = 205ELSEIF (IMG(NI,NJ).EQ.210) THENIMG(NI,NJ) = 220ENDIFENDIFN1N = SPH(NK,NL) + 1IF (IMG(NI,NJ).GT.0) THENNSPH(NIN) = NSPH(NIN) + 1ENDIF40 CONTINUEIF (NSPH(NSZ).LT.CHKSPH(NSZ)) GOTO 1120 CONTINUECC THIS IF BLOCK ENSURES THAT ALL POSSIBLE PORE SIZE VALUESC ARE CHECKED.C11 IF (L.GT.99) THENL= 101ELSEL=L+2ENDIFCC THIS LOOP DETERMINES THE MAXIMUM SIZE OF CiRCLE WHICH ISC COMPLETELY WITHIN THE PORE OR PORE THROAT.CDO 50 I=2,LIF (NSPH(I).EQ.CHKSPH(I)) THENRFIT = I-IENDIF50 CONTINUECC THE VALUE OF NREF IS DETERMINED: IF THE CENTER USED IN THISC SUBROUTINE WAS THE CENTER OF A PORE, NREF IS SET TO 255, AN])C IFITWASTHECENTEROFAPORETHROAT,ITISSETTO2O.CDO 60 I NI1-1,NIl+1DO 60 J=NJ1-l,NJ1+lIF (IMG(I,J).EQ.255) THENNREF = 2551 63ELSEIF (IMG(I,J).EQ.20) NREF =20ENDIF60 CONTINUECC THIS FINAL CHECK ENSURES THAT ANY CIRCLES WHICH TERMINATE AGAINSTC THE EDGES OF THE IMAGE ARE NOT INCLUDED IN THE PORE SIZEC DISTRIBUTION, AS THESE PORES ARE NOT COMPLETELY INCLUDED WITHiNC THE IMAGE.CIF ((NI1-RFIT).LT.4) THENNREF =0ELSEIF ((NI1+RFIT).GT.515) THENNREF=0ELSEIF ((NJ1-RFIT).LT.4) THENNREF=0ELSEIF ((NJ1+RFIT).GT.771) NREF =0ENDIFENDIFENDIFRETURNEND164APPENDIX 3Berea 100 ConfocalPore SizeDistributions1 65Berea 100 Pore Size Distributions0.120.080.040.001060.120.080.040.00Image ABR1ii’Image BBR1II0QE00010Pore Radius (m)Aio5Pore Radius (m)1 66Berea 100 Pore Size Distributions0.100.080.060.040.020.000.081060.060.040.020.00 -Image CBR1Image DBR110110000C.)E010Pore Radius (m)—•1Iio5Pore Radius (m)1 670I.C)0Image FBR1L10Berea 100 Pore Size Distributions0.100.080.060I.C.)I 0.040.0210-sPore Radius (m)Image GBR10.150.10-0.050.00 -I106LLi04Pore Radius (m)1 68Berea 100 Pore Size Distributions0.080.060.040.020.0010.6Image HBR1Image JBR1I I’CC)IC•—C)Pore Radius (m)0.100.080.060.040.020.00106Pore Radius (m)1 69Berea 100 Pore Size Distributions0.080.060.040.020.00 -11060.12 -0.08Image KBR1Image LBR110I0QJ0•—‘SJPore Radius (m)0.040.0010556IPore Radius (m)1 700C)E00.—IS011)0Berea 100 Pore Size DistributionsImage MBR1‘I’llrhPore Radius (m)0.080.060.04’0.020.001060.08 —0.060.040.020.00106Image NBR1Pore Radius (m)101710C)a)E0CC.)II jBerea 100 Pore Size DistributionsImage OBR10.100.080.060.040.02-0.00 - —1060.120.080.04Pore Radius (m)Image PBR10.00106‘ii...10Pore Radius (m)1 72Berea 100 Pore Size Distributions0.00-1060.080.060.040.020.00Image QBR1Image RBR10.120.080.04iIlA0.—a)00.—C)010Pore Radius (m)10I10Pore Radius (m)1 73APPENDIX 4S55 Confocal PoreSize Distributions1 74Carbonate (S55) Pore Size DistributionsImage A780.2 —CC.).) 0.1•00.0 -1060.120.—010Pore Radius (m)10Image B780.080.04’Pore Radius (m)1 750.-4‘-IC’SE00C.)C’SJCarbonate (S55) Pore Size DistributionsImage C780.1510Pore Radius (m)0.10•0.050.00•0.150.100.050.00•106Image D78Pore Radius (m)101 760.0Carbonate (S55) Pore Size DistributionsImage E780.15 -0.10-0.05 -0.00 41060.3LPore Radius (m)Image F780000.20.110Pore Radius (m)10177Carbonate (S55) Pore Size DistributionsICISC0.20.3Image G78Image H78100.110Pore Radius (m)0.20.1•IPore Radius (m)1 780.—-‘S0‘-4E0Carbonate (S55) Pore Size DistributionsImage J780.12 —00.080.040.00 -.1060.20.1’10Pore Radius (m)Image K780.0•Pore Radius (m)10’1 790C)1)E0I0.2Pore Radius (m)LCarbonate (S55) Pore Size DistributionsImage L780.2 —0.10.0106 10Pore Radius (m)Image M780.11 80Carbonate (S55) Pore Size DistributionsCCI0.2Image N78Image 078100.1-1Li05Pore Radius (m)CCQC0.20.1’0.0106 10Pore Radius (m)181CIjACarbonate (S55) Pore Size DistributionsImage P780.12CI.C.)j0.080.04io-5Pore Radius (m)10Image Q780.120.080.040.00106Pore Radius (m)1 82Carbonate (S55) Pore Size DistributionsImage R780.20.1•0.0•Image S7810C0I’0I0.21io4C000)E0io5Pore Radius (m)0.1•Pore Radius (m)1 830.—0.—C.)IACarbonate (S55) Pore Size DistributionsImage T780.10 —0.080.060.040.020.001060.30.20.1•0.0Pore Radius (m)i04Image U78106 10Pore Radius (m)io41 84APPENDIX 5S77 Confocal PoreSize Distributions1 85Carbonate (S77) Pore Size DistributionsImage Fl0.3 -:z0.0- II1106 io io io3Pore Radius (m)Image Gi0.50.40.3J 0.2106 io4 io3Pore Radius (m)1 86Carbonate (S77) Pore Size DistributionsImage HiIPore RadiusC0C.)IC0.—C)I1 •(m)0.40.30.20.1•1060.6 —0.50.40.30.20.1•0.010Image Ii106Pore Radiusio(m)101 870C.)j0.—C.)JCarbonate (S77) Pore Size DistributionsImage Ji0.40.30.20.1I0.0 —106I10Pore Radius (m)Image Ki0.30.20.10.0106 10Pore Radius (m)1 880.—00IImage Li106Carbonate (S77) Pore Size Distributions0.20.10.30.20.1•0.0•106io-4 io-3Pore Radius (m)Image MiIAi05 io4Pore Radius (m)101 890.,-I.io-5Pore Radius (m)Carbonate (S77) Pore Size DistributionsImage Ni0.3 -________0.2 -0.10.0 —‘1060.30.20.1•o•o• r.106Pore Radius (m)Image Oi0I 11101 90Carbonate (S77) Pore Size DistributionsC()ECImage P1io5 1o 1oPore Radius (m)0I.I0.20.1•1060.20.1•0.0• r106Image Q110Pore Radius (m)191Carbonate (S77) Pore Size Distributions0.—C)0Image Ri0.30.20.10.0• • rrr106 io i-Pore Radius (m)i031 92APPENDIX 6Carbonate SamplesNMR T1 and PoreSize Distributions193SAMPLE S31Continuous Relaxation Spectrum0.080.06{0.020.00 -1 10 100 1000Spin-Lattice Relaxation Times (ms)Discrete Relaxation Spectrum0.6050.4< 0.30.2-0.10.0 -1 10 100 1000Spin-Lattice Relaxation Times (ms)1 94SAMPLE S31Pore Size SpectraFast Diffusion Method0.080.060.040.020.00 •..108 10.6 ioPore Radius (m)Pore Size SpectrumSlow Diffusion Method0.3C-4t0.2A106 10_s ioPore Radius (m)1 95SAMPLE S32Continuous Relaxation Spectrum0.08 -0.06—0.040.020.00 I .....•1 10 100 1000Spin-Lattice Relaxation Times (ms)Discrete Relaxation Spectrum0.6050.40.30.2____,__I_,,FTTDTflI..........,. I1 10 100 1000Spin-Lattice Relaxation Times (ms)1 96SAMPLE S32Pore Size SpectraFast Diffusion Method0.08o 0.060.04C0.020.001O8 106 io4Pore Radius (m)Pore Size SpectrumSlow Diffusion Method0.2 -CI.Qa) 0.1C0.0 Ii0 106 io5 ioPore Radius (m)1 97SAMPLE S33Continuous Relaxation Spectrum0.08EA0.00 .. I1 10 100 1000Spin-Lattice Relaxation Times (ms)Discrete Relaxation Spectrum- 05j< 0.3a)0.2.... .1 10 100 1000Spin-Lattice Relaxation Times (ms)1 980.—I.C.)C0ElSAMPLE S33Pore Size SpectraFast Diffusion Method0.080.060.040.020.00iO8 iOPore Radius (m)Pore Size SpectrumSlow Diffusion Method0.30.20.10.0*.J Ii07 ;(.6Pore Radius (m)101 99SAMPLE S34Continuous Relaxation Spectrum0.08-1 10 100 1000Spin-Lattice Relaxation Times (ms)Discrete Relaxation Spectrum0.6a)05.0.4W0.30.20.1-0.0- —‘-—=1 10 100 1000Spin-Lattice Relaxation Times (ms)2000.—0000ISAMPLE S34Pore Size SpectraFast Diffusion Method0.080.060.040.020.000..8 iOPore Radius (m)Pore Size SpectrumSlow Diffusion Method0.30.20.10.0*j\.10 106Pore Radius (m)10201SAMPLE S35Continuous Relaxation Spectrum0.080.060.040.02o.00 I1 10 100 1000Spin-Lattice Relaxation Times (ms)Discrete Relaxation Spectrum0.6a)050.40.3wa)0.2W.=.. ...1 10 100 1000Spin-Lattice Relaxation Times (ms)2020C)0SAMPLE S35Pore Size SpectraFast Diffusion Method0.080.060.040.020.00i0..8 10.6 °Pore Radius (m)Pore Size SpectrumSlow Diffusion Method0C)J0.40.30.20.10.0*..Iio7 106Pore Radius (m)203SAMPLE S36Continuous Relaxation Spectrum0.080.060.001 10 100 1000Spin-Lattice Relaxation Times (ms)Discrete Relaxation Spectrum0.60•0.40.30.20.10.0 ==c—,—,—,-’-,-1 10 100 1000Spin-Lattice Relaxation Times (ms)2040.—c)C0I.c)1)ECSAMPLE S36Pore Size SpectraFast Diffusion Method0.08 -0.060.040.020.00i08 OPore Radius (m)Pore Size SpectrumSlow Diffusion Method0.30.20.10.0*..106 io4Pore Radius (m)205SAMPLE S37Continuous Relaxation Spectrum0.080.060.00•1 10 100 1000Spin-Lattice Relaxation Times (ms)Discrete Relaxation Spectrum0.6ti)050.4030.2.=1 10 100 1000Spin-Lattice Relaxation Times (ms)206SAMPLE S37Pore Size SpectraFast Diffusion Method0.08o 0.06O.0400.020.00-108 io 106 io5 ioPore Radius (m)Pore Size SpectrumSlow Diffusion Method0.300.20.1-. Ai07 106Pore Radius (m)207SAMPLE S38Continuous Relaxation Spectrum0.080.060.040.020.001 10 100 1000Spin-Lattice Relaxation Times (ms)Discrete Relaxation SpectrumSpin-Lattice Relaxation Times (ms)a)a).—a)0.60.50.40.30.20.10.0 GWQJ.. oooa1 10 100 1000208SAMPLE S38Pore Size SpectraFast Diffusion Method0.08i08 io7 i06 io io4Pore Radius (m)Pore Size SpectrumSlow Diffusion Method0.300.2I0.1•.A..1i07 106 io5 io4Pore Radius (m)209SAMPLE S39Continuous Relaxation Spectrum0.08I1 10 100 1000Spin-Lattice Relaxation Times (ms)Discrete Relaxation Spectrum0.6I)00.40.3a)0.20.1•0.0• DC1 10 100 1000Spin-Lattice Relaxation Times (ms)210SAMPLE S39Pore Size SpectraFast Diffusion Method0.080.060.040.020.000.—cj00I.1)$0108Pore Radius (m)Pore Size SpectrumSlow Diffusion Method0.30.20.10.0*../Pore Radius10(m)10211SAMPLE S40Continuous Relaxation Spectrum0.080.06 -0.040.020.001 10 100 1000Spin-Lattice Relaxation Times (ms)Discrete Relaxation SpectrumSpin-Lattice Relaxation Times (ms)Ia)a)..a)0.80.60.40.20.0 .DD1 10 100 10002120.08SAMPLE S40Pore Size SpectraFast Diffusion Method0.06.—C.)c0.04,C> 0.020.0010.8 O 10.6Pore Radius (m)Pore Size SpectrumSlow Diffusion Method0C.)110.30.20.10.0__106Pore Radius (m)213SAMPLE S41Continuous Relaxation SpectrumO.080.060.040.020.00 I1 10 100 1000Spin-Lattice Relaxation Times (ms)Discrete Relaxation Spectrum0.6a0.4< 0.3p. 0.20.1•0.0 ;I.fl.I.•.I.fl.1fl. i:i:i;ixR—————r——r——i—t—i—n1 10 100 1000Spin-Lattice Relaxation Times (ms)2140.08SAMPLE S41Pore Size SpectraFast Diffusion Method0.060.040.020.000I.c)00I’0108 io7 106 io ioPore Radius (m)Pore Size SpectrumSlow Diffusion Method0.30.20.10.0ioPore Radius10(m)10215SAMPLE S42Continuous Relaxation Spectrum0.080.060.040.020.00• I I1 10 100 1000Spin-Lattice Relaxation Times (ms)Discrete Relaxation Spectrum0.6a)a0.4< 0.3a)0.20.1•0.0 .IiI000aaDo1 10 100 1000Spin-Lattice Relaxation Times (ms)216SAMPLE S42Pore Size SpectraFast Diffusion Method0.080.060.040.020.00CI.C)1CCI0I108Pore Radius (m)Pore Size SpectrumSlow Diffusion Method0.30.20.10.0_/.Aiø ioPore Radius (m)217SAMPLE S43Continuous Relaxation Spectrum0.080.060.040.020.00-1 10 100 1000Spin-Lattice Relaxation Times (ms)Discrete Relaxation Spectrum0.60.4• 0.3JL1 10 100 1000Spin-Lattice Relaxation Times (ms)218SAMPLE S43Pore Size SpectraFast Diffusion Method0.080.060.040.020.00CCI.c)CC0QCi08 iOPore Radius (m)Pore Size SpectrumSlow Diffusion Method0.30.20.10.010A•1116Pore Radius (m)102190.06 -0.040.02 -0.00 -Spin-Lattice Relaxation Times (ms)SAMPLE S44Continuous Relaxation Spectrum0.081 10 100 1000Spin-Lattice Relaxation Times (ms)Discrete Relaxation SpectrumI1)0.60.50.40.30.20.10.0 Looaa j...Lq..1 10 100 100022000SAMPLE S44Pore Size SpectraFast Diffusion Method0.080.060.040.02 -0.0010.8 °Pore Radius (m)Pore Size SpectrumSlow Diffusion MethodC0I0.30.20.10.0*..I10 1O6Pore Radius (m)i04221SAMPLE S45Continuous Relaxation Spectrum0.081 10 100 1000Spin-Lattice Relaxation Times (ms)Discrete Relaxation Spectrum0.5 -0.4 -.1 10 100 1000Spin-Lattice Relaxation Times (ms)2220C.)C0C.)zC>SAMPLE S45Pore Size SpectraFast Diffusion Method0.080.060.040.020.00iO8 0 0Pore Radius (m)Pore Size SpectrumSlow Diffusion Method0.30.20.10.0*.i .io6 ioPore Radius (m)223SAMPLE S46Continuous Relaxation Spectrum0.080.060.040.020.0o I1 10 100 1000Spin-Lattice Relaxation Times (ms)Discrete Relaxation Spectrum0.60.50.41 10 100 1000Spin-Lattice Relaxation Times (ms)224SAMPLE S46Pore Size SpectraFast Diffusion Method0.080.060.040.020.00C.—CCIi0.8 1O6 0Pore Radius (m)Pore Size SpectrumSlow Diffusion Method0.30.20.10.01A!i-Pore Radius (m)225SAMPLE S47Continuous Relaxation Spectrum0.080.060.04 -0.02 -0.00 -1 10 100 1000Spin-Lattice Relaxation Times (ms)Discrete Relaxation SpectrumSpin-Lattice Relaxation Times (ms)E—ci)0.80.60.40.20.0L1 10 100 1000226SAMPLE S47Pore Size SpectraFast Diffusion Method0.08o 0.060.0400.020.00 I108 iO-7 10 io5Pore Radius (m)Pore Size SpectrumSlow Diffusion Method0.30: ::0.0• Aio-7 io-4Pore Radius (m)227SAMPLE S48Continuous Relaxation Spectrum0.080.060.040.02 -0.00 - ——I1 10 100 1000Spin-Lattice Relaxation Times (ms)Discrete Relaxation SpectrumSpin-Lattice Relaxation Times (ms)I0.80.60.40.20.0 JDODDD1 10 100 10002280C)00.r-)ISAMPLE S48Pore Size SpectraFast Diffusion Method0.080.060.040.020.0010.8 O. °Pore Radius (m)Pore Size SpectrumSlow Diffusion Method0.30.20.10.0*..\.i07 ioPore Radius (m)10229SAMPLE S49Continuous Relaxation Spectrum0.080.06 -0.040.020.001 10 100 1000Spin-Lattice Relaxation Times (ms)Discrete Relaxation SpectrumSpin-Lattice Relaxation Times (ms)a0.80.60.40.20.0*..E.1 10 100 10002300.—C.)00C.)jSAMPLE S49Pore Size SpectraFast Diffusion Method0.080.060.040.020.000.30.20.10.01-108 10 106 io5 io4Pore Radius (m)Pore Size SpectrumSlow Diffusion Method10.6Pore Radius (m)i04231SAMPLE S54Continuous Relaxation Spectrum0.070.00 I10 100 1000 10000Spin-Lattice Relaxation Times (ms)Discrete Relaxation SpectrumI)050.4W0.3W......1 10 100 1000 10000Spin-Lattice Relaxation Times (ms)232SAMPLE S54Pore Size SpectraFast Diffusion Method0.070.06•, 0.05C)0.04000E 0.030.020.01-i07 106 IIYPore Radius (m)Pore Size SpectrumSlow Diffusion Method0.80.6C)0.4C0.20.0106 io5 io4Pore Radius (m)233SAMPLE 855Continuous Relaxation Spectrum0.07 -0.00 I I10 100 1000 10000Spin-Lattice Relaxation Times (ms)Discrete Relaxation Spectrum0.60-e 05-I____ _____0.0 — .....__r_..r_r_r.......rITrTrrlIdIrrrrrrr •.I -1 10 100 1000 10000Spin-Lattice Relaxation Times (ms)234SAMPLE S55Pore Size SpectraFast Diffusion Method0.060.050.040.030.020.010.000C.)00C.)C.)z0>102 1111 101Pore Radius (m)Pore Size SpectrumSlow Diffusion Method0.30.20.10.0106Pore Radius (m)10235SAMPLE S77Continuous Relaxation Spectrum0.080.060.040.020.0010 100 1000 10000Spin-Lattice Relaxation Times (ms)Discrete Relaxation Spectrum0.80.60.41 10 100 1000 10000Spin-Lattice Relaxation Times (ms)236SAMPLE S77Pore Size SpectraFast Diffusion MethodO.080.06g 0.04C0.020.00i0 10.6Pore Radius (m)Pore Size SpectrumSlow Diffusion Method0.60.50.—0.4) 0.3I0.20.1•IPore Radius (m)237APPENDIX 7NMR Ti distributionsand M(O) vs. Swgraphs for PartiallySaturated carbonatesampies238SampleS31NormalizedTiSpectraatPartialSaturations0.08Saturation0.9840.8570.7740.7190.6530.060.040.020.00—0--110100Spin-LatticeRelaxationTime(ms)1000F’.)CA)CD220020001800M(0)160014001200SampleS31CorrelationofSaturationwithM(O)0.60.70.80.9Saturation1.0F’)ISampleS41NormalizedTiSpectraatPartialSaturations0.080.060.040.020.00Saturation--0.9830.9200.8320.7960.6730.560110100Spin-LatticeRelaxationTime(ms)10001)1300120011001000 900800700600SampleS41CorrelationofSaturationwithM(O)M(0)0.50.60.70.80.9Saturation1.0F\)I\)SampleS42NormalizedTiSpectraatPartialSaturations0.08Saturation0.9890.8720.7110.6110.5120.060.040.020.00—0---110100Spin-LatticeRelaxationTime(ms)1000I’) c)SampleS42CorrelationofSaturationwithM(O)M(O)1200 800400 0.40.60.8Saturation1.0r’3 -a ISampleS44NormalizedTiSpectraatPartialSaturations0.060.050.040.030.020.010.00-0----0-Saturation0.9140.8120.7110.6140.5090.324110100Spin-LatticeRelaxationTime(ms)1000C,’0)M(O)16001200 800400SampleS44CorrelationofSaturationwithM(O)0.3Saturation0.50.70.9SampleS54NormalizedTiSpectraatPartialSaturationsSaturations0.9950.920—G—0.886—.———0.7940.730—0-——0.6320.5570.4640.362-I-0.237ro —1Spin-LatticeRelaxationTime(ms)ii)0.080.060.040.020.0010100100010000I.’)-700060005000M(0)4000300020001000SaturationSampleS54CorrelationofSaturationwithM(O)0.20.40.60.81.0F’)- (.0Spin-LatticeRelaxationTime(ms)—0--.1.—--ISampleS55NormalizedTiSpectraatPartialSaturations0.07 0.060.050.040.030.020.010.00Saturation0.9840.8700.7780.6970.4650.4050.3080.1840.035110100100010000SampleS55CorrelationofSaturationwithM(O)600050004000M(O)300020001000 00.00.20.40.60.8Saturation1.0I\)01 0

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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>
                        
                    
IIIF logo 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-0052496/manifest

Comment

Related Items