Physics Reports (
)
–
Contents lists available at SciVerse ScienceDirect
Physics Reports journal homepage: www.elsevier.com/locate/physrep
Magnetic resonance imaging in laboratory petrophysical core analysis J. Mitchell a,b,∗ , T.C. Chandrasekera a , D.J. Holland a , L.F. Gladden a , E.J. Fordham b a
Department of Chemical Engineering and Biotechnology, University of Cambridge, Pembroke Street, Cambridge CB2 3RA, United Kingdom
b
Schlumberger Gould Research, High Cross, Madingley Road, Cambridge CB3 0EL, United Kingdom
article
info
Article history: Accepted 9 January 2013 Available online xxxx editor: M.L. Klein Keywords: MRI Petrophysics Core analysis Heterogeneity mapping Oil recovery Capillary pressure
abstract Magnetic resonance imaging (MRI) is a well-known technique in medical diagnosis and materials science. In the more specialized arena of laboratory-scale petrophysical rock core analysis, the role of MRI has undergone a substantial change in focus over the last three decades. Initially, alongside the continual drive to exploit higher magnetic field strengths in MRI applications for medicine and chemistry, the same trend was followed in core analysis. However, the spatial resolution achievable in heterogeneous porous media is inherently limited due to the magnetic susceptibility contrast between solid and fluid. As a result, imaging resolution at the length-scale of typical pore diameters is not practical and so MRI of core-plugs has often been viewed as an inappropriate use of expensive magnetic resonance facilities. Recently, there has been a paradigm shift in the use of MRI in laboratory-scale core analysis. The focus is now on acquiring data in the laboratory that are directly comparable to data obtained from magnetic resonance well-logging tools (i.e., a common physics of measurement). To maintain consistency with well-logging instrumentation, it is desirable to measure distributions of transverse (T2 ) relaxation time – the industry-standard metric in well-logging – at the laboratoryscale. These T2 distributions can be spatially resolved over the length of a core-plug. The use of low-field magnets in the laboratory environment is optimal for core analysis not only because the magnetic field strength is closer to that of well-logging tools, but also because the magnetic susceptibility contrast is minimized, allowing the acquisition of quantitative image voxel (or pixel) intensities that are directly scalable to liquid volume. Beyond simple determination of macroscopic rock heterogeneity, it is possible to utilize the spatial resolution for monitoring forced displacement of oil by water or chemical agents, determining capillary pressure curves, and estimating wettability. The history of MRI in petrophysics is reviewed and future directions considered, including advanced data processing techniques such as compressed sensing reconstruction and Bayesian inference analysis of under-sampled data. Although this review focuses on rock core analysis, the techniques described are applicable in a wider context to porous media in general, such as cements, soils, ceramics, and catalytic materials. © 2013 Elsevier B.V. All rights reserved.
Contents 1. 2.
Introduction............................................................................................................................................................................................. NMR in petrophysics: fundamentals .....................................................................................................................................................
2 7
∗ Corresponding author at: Department of Chemical Engineering and Biotechnology, University of Cambridge, Pembroke Street, Cambridge CB2 3RA, United Kingdom. Tel.: +44 0 1223 325426. E-mail addresses:
[email protected],
[email protected] (J. Mitchell). 0370-1573/$ – see front matter © 2013 Elsevier B.V. All rights reserved. doi:10.1016/j.physrep.2013.01.003
2
J. Mitchell et al. / Physics Reports (
2.1.
3.
4.
5.
6.
7.
)
–
Relaxation time measurements ................................................................................................................................................. 2.1.1. Longitudinal relaxation ............................................................................................................................................... 2.1.2. Transverse relaxation .................................................................................................................................................. 2.1.3. Pore size determination .............................................................................................................................................. 2.2. Diffusion measurements ............................................................................................................................................................ 2.3. Internal gradients ....................................................................................................................................................................... 2.4. Multi-dimensional correlations ................................................................................................................................................. 2.4.1. nD inversions ............................................................................................................................................................... 2.4.2. Relaxation time correlations....................................................................................................................................... 2.4.3. Diffusion-relaxation correlations ............................................................................................................................... MRI in petrophysics: fundamentals ....................................................................................................................................................... 3.1. Definitions and terminology ...................................................................................................................................................... 3.2. Basic imaging theory .................................................................................................................................................................. 3.2.1. Frequency encoding .................................................................................................................................................... 3.2.2. Phase encoding ............................................................................................................................................................ 3.2.3. Slice selection............................................................................................................................................................... High resolution MRI in core analysis ..................................................................................................................................................... 4.1. Heterogeneity mapping — rapid imaging techniques .............................................................................................................. 4.2. Heterogeneity mapping — quantitative imaging techniques .................................................................................................. 4.3. Fluid transport ............................................................................................................................................................................ 4.3.1. Monitoring imbibition and recovery processes ......................................................................................................... 4.3.2. Flow mapping .............................................................................................................................................................. 4.4. Capillary pressure measurements ............................................................................................................................................. Extension to very low field strengths .................................................................................................................................................... 5.1. Relaxation time mapping ........................................................................................................................................................... 5.2. Comparison of data..................................................................................................................................................................... 5.3. Monitoring forced displacement of oil ...................................................................................................................................... 5.4. Capillary pressure measurements ............................................................................................................................................. 5.5. Multi-dimensional relaxation and diffusion mapping ............................................................................................................. Advanced imaging techniques ............................................................................................................................................................... 6.1. Compressed sensing ................................................................................................................................................................... 6.2. Grain sizing: a Bayesian-MR approach...................................................................................................................................... Conclusions.............................................................................................................................................................................................. Acknowledgments .................................................................................................................................................................................. Appendix A. Nomenclature ............................................................................................................................................................... Appendix B. Abbreviations ................................................................................................................................................................ References................................................................................................................................................................................................
7 7 8 9 10 11 12 13 14 14 16 16 17 17 18 18 19 19 22 25 27 28 31 34 34 36 39 41 43 45 45 47 50 51 51 53 54
1. Introduction Understanding the structure-transport relationships that govern the behavior of fluids in porous media is a complicated topic of significant industrial and academic interest [1]. Porous media are ubiquitous in nature and industry, and are encountered in many commercial applications including chemical reactors, construction materials, and underground reservoirs containing oil, gas, or water. The recovery of oil from reservoirs is the motivation behind the experimental techniques discussed in this review. The properties of both the fluids (oil, water or brine, miscible or immiscible gas) present in the reservoir and the rock matrix (pore network, mineral composition, mechanical performance) need to be assessed in order to predict remaining oil reserves and potential recovery accurately. As the global price of oil rises, enhanced oil recovery (EOR) methods such as polymer or CO2 injection are becoming economical. These techniques bring their own research challenges: for example, the understanding of non-Newtonian fluid behavior in porous media. Ideally, all the properties of nascent and injected fluids, plus geophysical parameters of the rock, would be determined in situ in the reservoir using a suite of well-logging tools. However, important measurements of multiple fluid-phase capillary pressure and permeability are impractical with current wireline technology. As such, laboratory-scale studies of cored reservoir material are a critical component in oil exploration, production, and monitoring. Nuclear magnetic resonance (NMR) is recognized by the petrophysical community in the form of well-logging tools: a single-sided permanent magnet and radio frequency (rf) resonator are lowered into a well in order to interrogate the fluid properties in the formation near the bore. The experiments that can be conducted on such devices are limited to measurements of transverse relaxation time (T2 ) and self-diffusion coefficient (D ). Additionally, current NMR logging tools can ‘‘see’’ only a few centimeters into the formation and as such the results can be dominated by the effects of drilling fluid invasion into the volume of rock studied. Notwithstanding, the data are often sufficient to distinguish fluid-phases and as the logging tool is drawn up the well, fluid saturation and porosity (total fluid volume) are determined as a function of depth, thereby providing a vertical spatial profile of these quantities. These well-logs have a typical vertical spatial resolution of a few decimeters. Some logging tools can also provide a radial (depth) profile into the formation. Accordingly, well-logging is considered a magnetic resonance imaging (MRI) technique. In this article we are concerned with MRI
J. Mitchell et al. / Physics Reports (
)
–
3
Fig. 1. A wide range of rock types are studied using MRI in core analysis. A selection of quarried rock core-plugs are shown here, including limestones (LS) and sandstones (SS). Specific examples are (back row, left to right) Doulting LS, Birchover SS, Ketton LS; (middle row, left to right) Bentheimer SS, Marble LS, Ohio Blue SS; (front row, left to right) Lazenby Red SS, and a different Birchover SS. Other rocks commonly studied in the literature include Fountainbleu SS, Berea SS, and Indiana LS. The rocks listed here are obtained from quarried formations and as such are well-consolidated (ideal) samples. Cores obtained from oil field reservoirs may include poorly consolidated sands and shales that are difficult to handle in the laboratory.
applications in laboratory-scale core analysis, and as such do not provide a detailed discussion on well-logging methodology. Comprehensive reviews of well-logging are to be found in a special issue of Concepts in Magnetic Resonance (volume 13, issue 6) with concise overviews in [2,3]. At the laboratory-scale, NMR is a powerful technique that provides information on fluids confined in porous materials. Understanding of oil recovery processes at the laboratory-scale is critical for generating accurate models to predict fieldscale performance in the reservoir. A variety of disciplines are encompassed by the broad term ‘‘NMR’’, the predominant being spectroscopy, imaging, and time domain analysis (spin relaxation and diffusion). All three concepts can be combined in a single measurement, although the hardware requirements of the acquisition modalities can be very different depending on the application. In this review we focus on MRI experiments appropriate for laboratory-scale core analysis. To provide sensitivity to spatial variations in porosity and fluid saturation, it is necessary to encode other parameters in the images, including chemical shift (spectroscopy), relaxation time, or diffusion coefficient. The NMR experiments available in the laboratory are far more diverse than those feasible down-hole. Under favorable conditions, the NMR signal intensity (number of resonant nuclei) is linearly proportional to fluid volume over a wide range of core-plug saturation states, allowing porosity and fluid-phase saturation to be determined. A T2 relaxation time distribution can be acquired and this is a standard measurement in NMR core analysis, being the metric of choice down-hole. T2 is sensitive to details of the confining geometry (pore size) [4] and the chemical species present, thus providing fluid-phase discrimination between oil and brine in the rock. Other forms of fluid-phase discrimination available in the laboratory include longitudinal relaxation time (T1 ) and D ; both of these metrics are also sensitive to pore geometry to differing extents. NMR spectroscopy provides the most direct fluid-phase discrimination as it is sensitive only to the chemical species present. The basic NMR toolbox is described briefly in Section 2. A powerful supplement to all these laboratory-scale measurements is the possibility of incorporating spatial dimensions with MRI [5]; these spatially resolved measurements are the focus of this review. The basics of MRI are covered in Section 3. A plethora of imaging techniques are available, tailored to suit almost any combination of application and sample — even in core analysis, sample properties can vary widely. Examples of some outcrop rocks used as proxies for reservoir material are shown in Fig. 1. The imaging sequences most appropriate for capturing details of millimeter-scale heterogeneity in core-plugs are described in Section 4. The combination of MRI and NMR techniques allows the core analyst to study petrophysical phenomena beyond bulk fluid-phase saturation or pore size distributions. Areas of current interest are the application of MRI to in situ monitoring of forced displacement of oil (including EOR methods), rapid capillary pressure curve measurements, and wettability alteration. Case studies highlighting these areas of interest are given in Section 5. Recent developments in advanced data processing algorithms, offering the possibility of time-resolved studies and determination of rock properties such as grain size distributions, are explained in Section 6. A wide range of NMR magnets are now available commercially and may be tailored to the requirements of the user. Spectrometers are often named after the field strength of their associated magnet and are classified using three broad terms: low-field, intermediate-field, and high-field. However, the exact values of B0 that each of these terms covers is not strictly defined. Typically, low-field refers to any magnetic field strength where chemical spectroscopy is impractical and encompasses permanent magnet technology. These systems are sold commercially for quality control in industrial environments and are used to measure relaxation times or diffusion coefficients. Accordingly, we define low-field to mean
4
J. Mitchell et al. / Physics Reports (
)
–
Fig. 2. Summary of magnetic field strengths — primary applications and measurements. The ranges of B0 encompassed by the terms low-field (B0 < 1 T), intermediate-field (1 T < B0 < 3 T), and high-field (B0 > 3 T) are indicated. There are no strict definitions of these ranges, and we have chosen the boundaries to correspond with the most common experiments conducted in each of these regimes: relaxation and diffusion (low-field for industrial process and quality control), imaging (intermediate-field for clinical diagnosis), and spectroscopy (high-field for molecular structure determination). In practice, relaxation, diffusion, and imaging may be conducted at any field strength, while spectroscopy is limited by magnetic field homogeneity which is not sufficient, typically, in low-field permanent magnets, although examples of low-field spectroscopy have been published [6]. We define a special subsection of the low-field regime as ‘‘very low field’’, referring to the typical field strengths of well-logging tools (B0 ∼ 10 mT); this is also the field strength preferred for laboratory-scale core analysis, although throughout the literature core-plug imaging has been performed at much higher field strengths. Photographs (top) show typical magnet designs encountered at the different field strengths: (left) an Oxford Instruments Geospec2 (B0 = 50 mT) benchtop permanent magnet, (middle) a Magnex horizontal bore superconducting ‘‘small animal’’ imaging magnet (B0 = 2 T), and (right) a Bruker vertical bore superconducting magnet (B0 = 7 T).
B0 < 1 T. This is a very broad category as the special cases of Earth’s field NMR (B0 ∼ 10 µT) and well-logging tools (B0 ∼ 10 mT) fall into this definition of low-field. In core analysis, it is important to differentiate bench-top magnets designed to be comparable to well-logging tools (B0 ≈ 50 mT) and typical bench-top spectrometers operating at B0 ∼ 0.5 T. We therefore introduce the concept of very low-field NMR, where 1 mT < B0 < 100 mT, corresponding to the usable range of magnetic fields on well-logging tools. Intermediate-field NMR covers the range of field strengths at which medical imaging for clinical diagnosis is typically implemented. High resolution images can be obtained at field strengths as low as B0 ∼ 0.5 T, with the upper range for full-body MRI scanners currently around B0 ∼ 4 T, although higher field full-body systems have been constructed. Functional magnetic resonance imaging (fMRI) may be conducted at higher field strengths. For the purposes of core analysis, we consider intermediate-field to span 1 T < B0 < 3 T. High-field NMR is considered the domain of spectroscopy, either solid or liquid state; high resolution spectra are obtained typically at field strengths of B0 > 3 T. Significantly higher field strengths are used for structural determination of complicated molecules, such as proteins. Commercial magnets are available up to B0 ≈ 22 T, corresponding to a 1 H Larmor frequency of ν0 = 950 MHz. The definitions of magnetic field strength are summarized in Fig. 2. Unfortunately, laboratory-scale MRI did not prove to be the core analyst’s panacea that was envisioned in the early 1990s. The majority of MRI studies of rocks have been implemented on intermediate to high-field magnets (B0 > 1 T) derived directly from the spectrometers available commercially for clinical diagnosis and biomedicine, as reviewed in [7]. The so-called ‘‘small animal imaging’’ class of magnets have been favored as they readily accommodate both core-plug and temperature/pressure regulated sample enclosure. In medical imaging, there is a constant drive to develop high-field magnet systems: a higher magnetic field provides better signal-to-noise ratio (SNR) and hence enables higher imaging resolutions, and improves chemical shift contrast. To some extent, the same trend has been followed in MRI of rock cores, including attempts to image pore-scale structures within the rock [8]. However, due to the magnetic susceptibility contrast between the solid and fluid which reduces the effective resolution, it is impractical to use high-field magnets (B0 > 3 T) to image the contents of individual pores. Additionally, magnetic susceptibility induced ‘‘internal gradients’’ have a significant influence on measurements of the diffusion-sensitive T2 relaxation time [9,10]; the magnitude of the internal gradients scales with B0 . The presence of strong internal gradients ensures the MRI signal obtained with conventional medical imaging technology is not quantitative, as discussed in Section 2.3, and therefore cannot be rescaled to fluid volume. The lack of quantification in intermediate to high-field imaging has been a stumbling block for MRI in core analysis and resulted in a steady decline in the use of such magnets, as shown in Fig. 3.
J. Mitchell et al. / Physics Reports (
)
–
5
Fig. 3. Historical survey of magnetic field strengths B0 (legend) employed for MRI of core-plugs. Early in the 1990s there was a preference for using horizontal bore superconducting ‘‘small animal’’ imaging magnets, typically with B0 = 2 T (intermediate-field) as these systems were accessible in many biomedical laboratories. However, the general lack of quantification in images of rocks obtained on such magnets lead to a rapid decline in the use of medical MRI systems for core analysis. MRI of rock core-plugs at other field strengths has remained reasonably constant since the first images published in the late 1970s. Recently, there has been an increase in low-field MRI publications with B0 < 1 T. This survey is intended to provide an indication of trends in MRI for core analysis, rather than be a comprehensive study; all publications surveyed are included in the reference list at the end of this article.
The specter of internal gradients is associated primarily with sandstone formations containing paramagnetic species [11]. Carbonates – either limestone (predominantly calcite: calcium carbonate, CaCO3 ) or dolostone (predominantly dolomite: calcium magnesium carbonate, CaMg(CO3 )2 ) – rarely contain significant quantities of paramagnetic ions, so the fluid content may be determined quantitatively at intermediate magnetic field strengths. Historically, there has been a preference for studying sandstones because of their microstructural simplicity, yet carbonates, notably in the Middle East, contain significant hydrocarbon reserves. Due to the limit on spatial resolution in three-dimensional (3D) imaging and unreliable quantification, MRI of core-plugs has often been viewed as an inappropriate use of expensive NMR facilities [12]. Despite these reservations, continued investment in low-field NMR hardware and technique development over the last decade means low-field MRI methods can be used to determine a range of petrophysical properties, including porosity, wettability, and capillary pressure. Accordingly, low-field MRI is now gaining popularity in core analysis. Another general limitation in NMR for core analysis is the sample volume that can be accommodated in commercial magnets. It is common to study cylindrical core-plugs of volume less than 100 cm3 , although custom NMR systems for scanning long cores have been developed [13,14]. When studying short plugs, non-uniform saturation can arise from microscopic or macroscopic heterogeneities, or fluid redistribution. In dynamic studies, end effects caused by geometry, capillarity, or viscous instabilities in the fluid flow are often significant in plugs of sub-decimeter length. The addition of spatial resolution enables these end effects to be discounted in the analysis of saturation states, as discussed in Section 5.3. Therefore, the physical limitations of current commercial magnets ensures imaging is an essential commodity, rather than a luxury, in NMR core analysis. We digress here to note that NMR scientists have introduced vagaries into the terminology used to describe the engineering aspects of core analysis. A rock ‘‘core’’ describes a cylindrical sample of formation material recovered from the reservoir; these may be full-diameter cores, extracted during drilling, with typical diameters ranging from 4.5 cm (1.75 in.) to 13.3 cm (5.25 in.) and up to 120 m (∼400 feet) in length. Sidewall-cores, extracted horizontally from the well-bore, are typically less than 2.5 cm (1 in.) in diameter and up to 4.5 cm (1.75 in.) in length. Given the range of sizes in which cores may be extracted, some can be accommodated by commercial NMR magnets. However, in core analysis it is more common to extract a cylindrical ‘‘core-plug’’ from a full-diameter core. These plugs are typically either 3.8 cm (1.5 in.) or 2.5 cm (1 in.) in diameter with a maximum length of about 7 cm (2.8 in.), although non-standard sizes may be cut, especially when the plug is manufactured to suit the dimensions of a NMR probe. Composite cores are constructed by stacking multiple plugs end-to-end in a core holder. Due to the variability in sample size, the distinction between cores and plugs should be emphasized when describing MRI experiments. Unfortunately, in some of the existing literature, the concept of a ‘‘core’’ has been applied loosely to mean any rock sample. The popularity of various MRI techniques has changed over the last three decades, driven by advances in experimental (hardware) design, understanding of the applicability of the methods, and general interest in different properties of rocks and fluids. Briefly, the main applications of MRI in core analysis are summarized as follows:
• Monitoring capillary adsorption. Low resolution one-dimensional (1D) profiling was introduced for monitoring capillary adsorption of water with examples of reservoir rocks (Locharbriggs sandstone, Corsehill sandstone, and Portland limestone) [15]. These measurements demonstrated the applicability of MRI to core-plugs and provided hydraulic diffusivities from capillary adsorption profiles. • Clay swelling and mud infiltration. Initial imaging studies conducted at intermediate-fields relied primarily on T1 mapping [16–19]; the T1 relaxation time is not sensitive to internal gradients [20]. These works included the application of MRI to the study of clays [16,17], which are of particular relevance to drilling mud formulation and mud invasion into
6
J. Mitchell et al. / Physics Reports (
)
–
the well bore, and in the study of shales [21]. Clay-bound water is notoriously difficult to measure directly because of the restricted molecular mobility. Dvinskikh et al. employed T1 contrast to enhance the sensitivity of the MRI measurement when monitoring swelling of bentonite clay [22]. Fagan et al. opted for continuous wave magnetic resonance imaging (CW-MRI) in order to probe the water in bentonite clay [23]. However, these specialist measurements have not been used widely in core analysis, although they offer insights into MRI techniques suitable for measuring rock samples with very high clay content. Details on relaxation time measurements are given in Section 2.
• Structure and flow mapping. Rapid (qualitative) imaging techniques developed for medical MRI have been used to determine structure and fluid transport in core-plugs: specifically, π -echo planar imaging (PEPI) [24,25], and later, rapid acquisition with relaxation enhancement (RARE) [26]. These high resolution 3D images are used to investigate the presence of sub-millimeter macroscopic structural heterogeneities in core-plugs, such as fossils or vugular porosity, that cannot be observed readily by other non-destructive means. Details of rapid heterogeneity mapping methods are given in Section 4.1. Further, these rapid imaging methods can be used to determine displacement maps in flowing systems, as will be discussed in Section 4.3.2.
• Porosity and saturation mapping. Quantitative images were first obtained through the acquisition of T2 relaxation time maps [27,28] and provide an accurate measure of local porosity (total signal intensity) and fluid saturation. Details on quantitative imaging techniques are given in Section 4.2. The addition of the T2 dimension provides sensitivity to fluidphases, making the experiment ideal for monitoring forced displacement of oil, see Section 4.3. Wettability can also be inferred from relaxation time measurements [29] as will be demonstrated in Section 5.5.
• Capillary pressure curve measurement. Quantitative intermediate-field imaging has been achieved using the single point imaging (SPI) [30–32] suite of methods to overcome the relaxation weighting in rapid medical imaging techniques, see Section 4.2. These techniques are appropriate for saturation profiling and have been applied to measurements of capillary pressure curves [33–37], building on previous work using other imaging methods [18,38–41]. Details on capillary pressure measurements by MRI are given in Section 4.4. There is current interest in extending the SPI class of measurements to probe local permeability [42,43] and provide a rapid determination of relative permeability. An alternative method of determining local permeability is gas imaging as the NMR properties of gases (spin density, relaxation time) depend on pressure [44]. MRI has also found niche applications in the oil and gas industry. For example, Coussot et al. used rheo-NMR to determine the complicated rheology in drilling muds [45]. For this application, a cone-and-plate rheometer (or Couette cell) was mounted in the magnet and the velocities of the rotating liquid measured as a function of shear rate. A review of rheoNMR is given in [46]. More unusual applications of MRI in core analysis include monitoring foam transport [47,48] and determining the oil saturation in fluid-phase sensitive pressure taps [49]. Recently, there has been a paradigm shift in the use of MRI in core analysis. Very low magnetic fields (B0 < 100 mT) are now preferred in order to provide comparability with well-logging tools. The use of low-field magnets also has the advantage of reducing the influence of internal gradients on measurements of T2 , allowing a wider range of rock types to be studied. The shift to very low-field experiments is accompanied by an inherent reduction in SNR compared to intermediatefield spectrometers. The SNR achievable at lower magnetic field strengths ensures that the maximum imaging resolution in heterogeneous materials is on the order of 100 µm (sub-millimeter) [50]. However, it is not practical, especially on low porosity rocks (φ ≤ 0.05), to divide the sample into voxels on the order of a cubic millimeter and still expect to achieve a usable SNR within a reasonable time. It is therefore more useful to focus on the acquisition of fast and quantitative 1D profiles. For example, T2 mapping on a low-field bench top spectrometer provides data that are directly analogous with a well-log: the variation of T2 with position (height) in the formation. The main difference between the laboratory-scale and field-scale measurements is in the spatial resolution. At the laboratory-scale, the typical resolution is on the order of millimeters, whereas the well-log resolution is on the order of decimeters. Imaging at low-field is not a new idea. The first full-body MRI scanners operated at low magnetic field strengths; for example, B0 = 0.04 T (1.7 MHz for 1 H) [51]. Low-field bench-top imaging magnets (B0 ∼ 0.5 T, νo ∼ 20 MHz) are currently popular in the pharmaceuticals industry for in situ monitoring of drug tablet swelling and assessment of coating efficacy [52]. Recently, a very-low field portable magnet (B0 ≈ 20 mT, ν0 = 1 MHz) was constructed for 2D imaging of trees to determine water transport through the trunk [53]. Many other mobile and low-field magnet designs exist, some of which are tailored for high-resolution profiling applications [54,55]. In Section 5 we focus on the state-of-the-art in very low-field MRI as applicable to core analysis. A 1D profile is sufficient for determining the saturation of plugs during forced drainage or imbibition because fluid transport occurs predominantly along one axis [56]. Spatially-resolved T2 distributions provide fluid-phase discrimination when monitoring forced displacement of oil by brine. When T2 contrast is insufficient for fluid-phase discrimination, two-dimensional (2D) D –T2 ‘‘diffusion editing’’ correlations [57] or T1 –T2 relaxation correlations [58] may be employed instead. These 2D acquisitions can be spatially resolved to provide additional insights [59]; for example, spatially-resolved T1 –T2 correlations provide an indication of wettability variation along a core-plug and are therefore useful when studying oil recovery by wettability altering agents. T2 maps can be used to infer pore body-to-throat ratios (BTR) when combined with MRI-centrifuge capillary pressure curve measurements. Overall, we demonstrate that MRI is implemented readily in very low-field experiments and offers significant advantages over bulk NMR analysis of core-plugs.
J. Mitchell et al. / Physics Reports (
)
–
7
To demonstrate the capabilities of very low-field instrumentation in core analysis applications, we present a series of experimental case studies throughout Section 5. These case studies consist of:
• Comparison of T2 maps obtained from doped water imaging phantoms using different imaging techniques. The results obtained with the different relaxation time mapping protocols available at low-field are compared in Section 5.2.
• Forced displacement of oil by brine is demonstrated in Section 5.3. The results for a model system of dodecane and brine in sandstone are described in detail, followed by a review of low-field oil recovery monitoring studies, including an example of chemical EOR. • The measurement of an air/brine capillary pressure curve for a limestone plug is demonstrated in Section 5.4. The chosen limestone has a bimodal pore size distribution, and the ability to select capillary pressure curves associated with either the microporosity or macroporosity through the utility of the T2 dimension is discussed. • An example of a spatially resolved, 2D relaxation time (T1 –T2 ) correlation is presented in Section 5.5. A composite plug of oil wet and water wet limestone is used to highlight the application to wettability mapping. Additionally, spatially resolved D-T2 correlations are reviewed. Finally, we consider very recent applications of advanced image processing techniques to core analysis. In this review we consider compressed sensing image reconstruction as a method of reducing acquisition times or improving resolution [60,61], see Section 6.1. Improvements to image acquisition speed enable more accurate studies of dynamic systems, including imbibition, drainage, and oil recovery. Another novel technique reviewed here is Bayesian inference analysis of MR data to determine grain size distributions [62], see Section 6.2. The ability to determine pore size, grain size, and flow on the same sample by combining several MRI techniques will provide important insights into the relationships between rock structure and fluid transport. 2. NMR in petrophysics: fundamentals To understand the techniques employed in MRI, it is necessary to have some knowledge of the basic principles of NMR. Magnetic resonance relies on the precession of nuclear spins when subjected to a strong magnetic field. The precession frequency, or Larmor frequency, ω0 corresponds to the radio frequency (rf) range of the electromagnetic spectrum and depends on the magnetic field strength as ω0 = γ B0 , where γ is the nucleus-specific gyromagnetic ratio. Following the usual convention, the direction of the static magnetic field is along the z-axis. By exciting the spins away from their equilibrium state on the z-axis with the application of a 90° rf pulse (equivalent to an applied magnetic field B1 ), it is possible to detect the spin precession in the x–y plane using a resonator tuned to the Larmor frequency. The applied rf is usually quoted as ν0 [MHz], where ν0 = ω0 /2π . The majority of NMR and MRI measurements on liquids in porous rocks focus on acquiring signal from 1 H nuclei. NMR is an insensitive measurement (detecting only a small net polarization) and so 1 H nuclei are advantageous in being approximately 100% naturally abundant and having a large γ = 2.65 × 108 rad T−1 s−1 . Accordingly, 1 H are the only nuclei detectable with current down-hole logging tool technology. Nevertheless, other nuclei are studied under laboratory conditions, typically with intermediate to high-field magnets where SNR is less of a concern. Carbon-13 offers the advantage of a much wider range of chemical shifts (relative to 1 H) and so is useful for identifying oil components. Additionally, 13 C allows oil to be detected in the presence of water when 1 H spectroscopy does not provide sufficient sensitivity [63], and may enable studies of hydrogen-free species such as CO2 . However, the spin-half carbon isotope has low natural abundance (∼1.1%) and a small γ = 6.73 × 107 rad T−1 s−1 so detection is difficult. The SNR in 13 C experiments may be improved through isotopic enrichment or polarization transfer techniques [64]. The other nucleus of importance in core analysis is 23 Na, useful for studying high salinity brines [65,66]. Despite the alternatives, 1 H remains the nucleus of choice in core analysis for the advantages noted already; therefore, unless otherwise stated, all experiments in this review pertain to 1 H NMR. 2.1. Relaxation time measurements Relaxation times are used to characterize pore size, wetting state, and fluid content in core-plugs. We review the basic techniques for determining longitudinal T1 and transverse T2 relaxation times. The numerical inversion required to convert relaxation time data to pore size distributions is discussed in Section 2.1.3. These relaxation time measurements are often combined with MRI techniques to provide relaxation weighted images or relaxation time maps. 2.1.1. Longitudinal relaxation Once excited, the spin ensemble will return to its equilibrium state on the z-axis following an exponential recovery. The time constant of this longitudinal recovery process is T1 . Accordingly, the recovery of signal (magnetization) b at time t is described by b (t ) beq (∞)
t = 1 − c exp − , T1
(1)
where c is a scalar determined by the starting condition of the spin ensemble and beq is the equilibrium magnetization corresponding to the fully recovered spin ensemble observed as t → ∞. As the NMR experiment is sensitive only to
8
J. Mitchell et al. / Physics Reports (
)
–
Fig. 4. Pulse sequence schematic for T1 relaxation time measurements via (a) inversion recovery and (b) saturation recovery [67]. In all cases, the thin and thick vertical bars represent 90° and 180° rf pulses, respectively. The saturation comb of 90° pulses is repeated nc times. Homospoil gradients (gray rectangles) may be included between the rf pulses to eliminate unwanted coherence pathways. The initial magnetization b at time t = 0 is set to (a) b(0) = −beq or (b) b(0) = 0. The spin ensemble is allowed to recover on the z-axis for a time τ1 before being interrogated by an excitation pulse. The observed FID is Fourier transformed to obtain a chemical spectrum. Alternatively, when chemical resolution is not required, the signal amplitude at time t = τ1 is determined simply from the initial points in the FID. The experiment is repeated for a range of τ1 times to determine the T1 recovery. A rapid double-shot T1 measurement is described in [68]; at present, this pulse sequence has not been used in combination with spatial resolution nor is it appropriate for very low-field applications due to the reliance on small tip angle rf pulse (reduced SNR).
magnetization in the x–y plane, it is necessary to excite the recovered spins away from the z-axis in order to determine the rate of recovery. The standard T1 measurement is described by Vold et al. [67]. The recovery to equilibrium is typically monitored by starting with the spin ensemble aligned along −z (inversion recovery, where c = 2) or saturated on the x–y plane (saturation recovery, where c = 1). Pulse sequence schematics for both methods are shown in Fig. 4. Inversion recovery provides twice the dynamic range in the data compared to saturation recovery, but relies on the application of a good 180° inversion rf pulse. Furthermore, the spin ensemble must be aligned along +z prior to the initial inversion pulse. To guarantee this condition, the recycle delay between scans tRD must be long and satisfy tRD ≡ 5 × T1 . Saturation of the spin ensemble is easier to achieve when B0 and B1 are non-uniform. Additionally, the saturation recovery experiment does not require a long tRD as saturation is achieved regardless of the initial state of the spin ensemble. The saturation recovery experiment can therefore be implemented much faster than the inversion recovery experiment, especially when long T1 times are observed (e.g., when T1 ∼ 10 s). 2.1.2. Transverse relaxation The observed NMR signal amplitude decays over time as the spin ensemble loses phase coherence in the x–y plane due to local magnetic field fluctuations. These fluctuations arise from dipolar interactions, heterogeneities in the background magnetic field, and other terms in the nuclear spin Hamiltonian. The observed signal decay, or free induction decay (FID) [69], is governed by the exponential time constant T2∗ , as [5] 1 T2∗
≈
1 T2
+ γ ∆B0 ,
(2)
where ∆B0 is a characteristic measure of magnetic field variation. T2∗ is inversely related to the full-width half-maximum (FWHM) of the spectral line, such that 1 T2
+ γ ∆B0 ≈ π FWHM.
(3)
When using superconducting magnets (typical in intermediate to high-field experiments), ∆B0 ≈ 0 in bulk liquids and so T2∗ ≈ T2 , where T2 is the transverse relaxation time determined predominantly by dipolar coupling, and J-coupling where relevant (e.g., in oils). However, when using permanent magnets, known for typically poor homogeneity compared to superconducting magnets, ∆B0 is large and the signal decay will be characterized by T2∗ ≪ T2 . A similar situation arises when measuring heterogeneous samples — such as liquids confined in porous media — even when B0 is homogeneous, because the magnetic susceptibility difference between the solid and liquid distorts the static magnetic field. This distortion gives rise to ‘‘internal gradients’’ within the pores; further details are given in Section 2.3. To measure the T2 relaxation in the presence of magnetic field inhomogeneities, the spin ensemble is refocused on the x–y plane to form a spin echo [70] by a 180° rf pulse applied a time τ2 after the initial excitation pulse. The spin echo, comprising two FIDs back-to-back, is centered on time te = 2τ2 after the excitation pulse. By repeatedly applying 180° rf pulses, separated in time by te , a train of echoes is formed [71]. The method of forming this echo train has been refined into the standard Carr–Purcell–Meiboom–Gill (CPMG) experiment [72]. The envelope of the echo maxima decays exponentially with T2 as b (t ) b (0)
= exp −
t T2
.
(4)
This pulse sequence allows the transverse relaxation time to be measured even in the presence of an inhomogeneous magnetic field, although errors can be introduced in the T2 measurement by off-resonance spin dynamics [73]. A schematic
J. Mitchell et al. / Physics Reports (
)
–
9
Fig. 5. Schematic of the standard CPMG pulse sequence [71,72] for T2 relaxation time measurements. In all cases, the thin and thick vertical bars represent 90° and 180° rf pulses, respectively. The echo time te = 2τ2 ; the echo maximum forms a time τ2 after the 180° pulse. The repeated application of 180° refocusing rf pulses generates a train of n echoes. After n echoes, a FID (latter half of final echo) is acquired and Fourier transformed to provide a chemical spectrum. The process is repeated for a range of n. When chemical sensitivity is not required (or not possible, especially at low-field), only a few data are acquired at the top of each echo to determine the maximum signal intensity. Following this protocol, an entire echo train is interrogated in a single scan.
of the CPMG pulse sequence is shown in Fig. 5. Refinements to the experimental timing have been suggested for use in the grossly inhomogeneous magnetic fields of unilateral devices [74], such as well-logging tools, although such modifications are not generally required for commercial laboratory instruments. In homogeneous magnetic fields, the Fourier transform (FT) of the FID provides the chemical spectrum of the fluid being measured. Chemical shielding (J-coupling) moderates homonuclear dipolar coupling and so spins in different chemical environments precess at slightly different frequencies. The 1 H chemical spectrum will contain spectral lines associated with each 1 H environment (e.g. CH2 , CH3 , OH, etc.) thereby allowing identification of the fluid components if multiple species are present, such as oil and water. Many excellent texts are available on the subject of chemical spectroscopy; see for example [6,75,76]. In non-uniform magnetic fields, 1 H spectra are typically dominated by the line broadening caused by the ∆B0 term in Eq. (2) and so spectroscopy of large samples, such as liquid-saturated rocks, is not practical with current permanent magnet technology. Permanent magnets capable of spectral resolution over very small, bulk liquid samples (diameter ∼5 mm) are reviewed in [6]. 2.1.3. Pore size determination The relaxation time constants T1 and T2 are determined by a number of factors including the mobility and physical/chemical environment of the spins; T2 is also sensitive to molecular diffusion in the presence of magnetic field gradients. When fluids are confined in porous materials, the spins associated with molecules colliding with the pore walls undergo enhanced relaxation. This enhanced relaxation is triggered by a change in molecular mobility [77] and/or dipolar coupling with paramagnetic species and other active sites on the solid surface [11,78]. If the spins at the pore surface exchange with those in the bulk liquid on the time-scale of the NMR experiment [79] (τ1 for T1 and 2τ2 for T2 ) then the observed relaxation rate is modified according to the surface-to-volume ratio S /V of the pores [4,80,81]: 1 T1,2
=
1 (bulk)
T1,2
+ ρ1,2
Sp Vp
,
(5) (bulk)
assuming fast diffusion such that ρ1,2 Vp /D Sp ≪ 1 [79], where T1,2 is the bulk liquid T1 or T2 relaxation time, respectively, ρ1,2 is a surface relaxivity term associated with T1 or T2 processes, respectively, and Sp /Vp is the surface-to-volume ratio of the pore network. A characteristic pore size ℓs is defined, considered equivalent to pore radius; for spherical geometry, 3/ℓs = Sp /Vp . Determination of ρ1,2 is not straight-forward, so re-scaling of relaxation time to pore size is not readily quantitative. In many limestones, ρ1,2 is calibrated by determining the surface area from krypton (gas) adsorption measurements [82]; however, such an approach is not possible in sandstones as clays adversely influence the adsorption process, so other methods of calibration must be sought. Alternative NMR methods for determining pore size in rocks include cryoporometry and restricted diffusion (diffusometry). NMR cryoporometry relies on the depression in freezing/melting point of a confined liquid, relative to the bulk liquid, to provide a pore size distribution [83]. Spatially resolved cryoporometry measurements on core-plugs were demonstrated in [84], although the method is only suitable for probing pore sizes ℓs ≪ 1 µm (microporosity) so the general application in core analysis is limited. NMR cryoporometry is reviewed in [85]. Restricted diffusion measurements are based on the modification to the self-diffusion coefficient of a liquid when the mean diffusion path-length is limited by the confining geometry. Theoretically, restricted diffusion will probe pore sizes in the range 1 µm ≤ ℓs ≤ 10 µm [86–88], although the technique has been shown to be unreliable in real reservoir rocks due to the interconnectivity between pores [89]. Relaxation time is one of the few parameters that is sensitive to a wide range of pores sizes, although this sensitivity relies on an appropriate degree of surface relaxation. If ρ1,2 is small, as is often the case in carbonate formations [82], then the term ρ1,2 Sp /Vp becomes vanishingly small in macropores (ℓs > 10 µm). Under such circumstances (bulk)
T1,2 = T1,2 , so only the bulk liquid relaxation time is observed. To generate a relaxation time distribution (considered equivalent to a pore size distribution), it is necessary to invert the relaxation time data. The NMR data are described generally by the 1st kind Fredholm integral equation b(t ) b(0)
∞
K(u, t )f (u)d (log10 u) + ε(t ),
= −∞
(6)
10
J. Mitchell et al. / Physics Reports (
)
–
where b represents the acquired NMR data, t is the experimental timing, u is the distribution of fitting parameters (relaxation time or diffusion coefficient), K is a kernel function describing the expected form of the data: for NMR data in general, K = exp {−t /u}; ε describes the noise. Solving for f (u) is a difficult inverse problem, even on the log10 scale. Therefore, inversion of NMR data is achieved using regularization techniques [90]. The problem can be expressed in the vector–matrix form b = Kf, where K is a kernel matrix corresponding to the theoretical attenuation value for a given measurement time, b is the acquired data vector and f is the probability distribution (vector) that we wish to determine. A solution ˆf is obtained that satisfies
ˆf = arg min ∥b − Kf∥22 ,
(7)
f ≥0
where ∥. . .∥2 is the ℓ2 (Euclidean) norm of the vector. However, as K is ill-conditioned, unphysical solutions may be obtained. A stable, physical solution is obtained in the presence of noise by constraining ˆf to be non-negative, bounded by upper and lower limits of u, and smooth. Tikhonov regularization [91] uses a penalty function to bias against solutions of ˆf containing unphysical fluctuations. Hence, we solve
ˆf = arg min ∥b − Kf∥22 + λ ∥Lf∥22 ,
(8)
f ≥0
where ∥b − Kf∥22 is the residual and indicates how close f is to being a true solution to the physical problem. The second term in Eq. (8), λ ∥Lf∥22 , is the penalty function, and this controls the degree of smoothing applied to the solution; L is the operator representing the choice of smoothness criterion, such as the 2nd derivative. The degree of smoothing is controlled by λ, the regularization parameter, which may be chosen by a variety of methods [92]; a popular and robust technique is generalized cross validation (GCV) [93]. 2.2. Diffusion measurements Pulsed field gradient (PFG) diffusion measurements were introduced by Stejskal and Tanner [94,95]. A magnetic field gradient of amplitude g is applied for a time tδ to provide the spins with a phase shift proportional to their position r(0). After an observation time t∆ , a second gradient of equal duration but amplitude −g is applied to remove the phase shift. If the spins have not moved, their net phase shift is zero. However, if the spins have moved to a new position r(t∆ ), then the spins have a net phase shift equal to
ϕ = γ g [r(t∆ ) − r(0)] .
(9)
This phase shift can be simplified to ϕ = 2π q · R where q = γ tδ g/2π is the magnetization wave vector and R = r(t∆ )− r(0) is the displacement vector. The signal decay due to diffusion is then b(g )/b(0) = exp −4π 2 q2 D t∆ ,
(10)
where q is the magnitude of q and D is the diffusion coefficient. The diffusion coefficient is obtained by repeating the measurement for a range of g. If all the timings remain constant, the only contribution to the normalized signal decay is diffusion. The basic principles of PFG diffusion measurements are the same regardless of the acquisition method. However, the details of Eq. (10) will depend on the exact form of the experiment. For example, the pulsed gradient spin echo (PGSE) experiment [94] under the narrow-pulse approximation (tδ ≪ t∆ ) with rectangular gradient pulses [5] gives b(g )/b(0) = exp −γ 2 tδ2 g 2 D t∆ .
(11)
The basic PSGE pulse sequence is shown in Fig. 6(a). The maximum amount of signal obtained in the PGSE sequence will depend on the ratio of T2 to t∆ . In rocks, T1 > T2 (even at low magnetic field strengths [96]) and so it is desirable to use a pulsed gradient stimulated echo (PGSTE) experiment instead [95]. A stimulated echo [70] stores the spin ensemble on the z-axis during the observation time, so the signal decays with T1 , enabling longer observation times to be employed. The basic PGSTE sequence is shown in Fig. 6(b); Eq. (11) still applies to this pulse sequence. As well as the applied gradient pulses, internal gradients within the pores will contribute to diffusive attenuation. In the presence of a background gradient g0 , the signal attenuation observed in a PGSTE experiment will be [97] b(g ) b(0)
=
1 2
2
exp −γ D
2
tδ
t∆ −
tδ 3
g2
2t 2 τse 2 + tδ 2τse t∆ − δ − tδ (tδ1 + tδ2 ) − tδ21 + tδ22 gg0 + τse2 t∆ − g0 . 3
3
(12)
Cotts et al. [97] devised several alternating pulsed gradient stimulated echo (APGSTE) sequences to reduce the influence of background gradients in diffusion experiments [97]. The most popular version is the Condition I ‘‘13-interval’’ APGSTE pulse
J. Mitchell et al. / Physics Reports (
)
–
11
Fig. 6. Basic diffusion pulse sequence schematics for (a) PGSE, (b) PGSTE, and (c) APGSTE. Thin and thick vertical lines indicate 90° and 180° rf pulses, respectively. Timings are specified following the notation of Tanner [95]: in each case, the spin echo time is τse and the observation time (defined as the time between leading gradient edges) is t∆ . Rectangular gradient pulses are shown with duration tδ and amplitude g. The gradient amplitude is incremented in successive experiments. Homospoil gradients are indicated by gray rectangles. The delay between the rf pulse and the leading gradient edge is tδ 1 , and the delay between the trailing gradient edge and the rf pulse is tδ 2 .
sequence, see Fig. 6(c). Generally, this 13-interval APGSTE experiment provides a signal b(g ) b(0)
=
1 2
exp −γ 2 D
tδ2 4
4t∆ − 2τse −
tδ 3
g 2 + τse tδ (tδ 1 − tδ 2 ) gg0 +
4 3
τse3 g02
,
(13)
such that the cross term in gg0 vanishes when tδ 1 = tδ 2 . Assuming g 2 ≫ g02 , the signal attenuation is dominated by the applied gradient pulses and not the internal gradients. The factor 12 in Eqs. (12) and (13) denotes that the stimulated echo returns half the initial signal amplitude. The data quality can be improved by using shaped gradient pulses (e.g., trapezoids or half-sine pulses) and such modifications to the experiment design result in an altered diffusion exponent. A summary of PFG experimental modifications and the associated expressions for signal attenuation are found in [98,99]. Improved compensation for background gradients can be achieved using the ‘‘17-interval’’ APGSTE pulse sequence [97] – or higher order versions [100] – which removes the cross term in gg0 regardless of the ratio of tδ 1 to tδ 2 and provides further reduction of the term in g02 relative to that in g 2 . However, the increased complexity of implementation and greater relaxation attenuation associated with these sequences means they are not widely used in core analysis. The pulse sequences described so far rely on fixed experimental timings and variable gradient amplitudes. However, in some well-logging tools the inherent magnetic field gradient in B0 resulting from the design of the unilateral magnet may be used as a diffusion gradient. The signal attenuation is achieved by increasing the spin echo time τse in a spin echo experiment, such that both t∆ and tδ are increased. In order to extract a diffusion coefficient, it is necessary to account for the increased relaxation time by deconvolving the contributions of D and T2 [57,101,102]. Further details on this experiment are provided in Section 2.4. An equivalent diffusion measurement, suitable for use in laboratory core analysis, has been published by Mitchell and Fordham [103]. The diffusion pulse sequences in Fig. 6 are also used to acquire a probability distribution of advective displacement (the so-called ‘‘flow propagator’’) [104], averaged over all possible starting positions of the spins in the volume of interest. The only modification required to the aforementioned pulse sequences is that the gradient amplitude is ramped between ±g in successive scans. Flow propagators have been used to assess pre-asymptotic Stokes’ flow in rocks [105–120]. Slice selection has also been combined with APGSTE flow propagator acquisitions to remove the influence of entry and exit effects from the sample on the displacement measurement [121]. 2.3. Internal gradients Magnetic susceptibility contrast induced internal gradients have a significant influence on measurements of T2 . The effective relaxation time T2,eff observed in the presence of a strong, spatially varying internal gradient will appear to contain multiple components due to the contributions from surface relaxation and diffusion in the CPMG data [10,122]. The internal gradients have no effect on T1 relaxation time measurements [20] and so the T1 distribution shape can be used as a comparator to assess the significance of diffusion on T2,eff . The strength of internal gradients scales with B0 but inversely with pore size ℓs . As spins diffuse through the internal gradients, enhanced dephasing of the spin ensemble √ occurs and this contributes to the loss in signal in a CPMG echo train. In bulk liquids, spins diffuse a distance ℓe ∼ D0 te , where D0 is the bulk (free) liquid diffusion coefficient. If the echo time te is sufficiently short, then the CPMG experiment compensates for diffusion. However, it is not practical to set te short enough in intermediate to high-field measurements of porous media to remove the sensitivity to diffusion in the internal gradients unless the effective gradient strength geff is
12
J. Mitchell et al. / Physics Reports (
)
–
Fig. 7. Comparison of porosity for a selection of sandstone core-plugs obtained via buoyant weight measurement and NMR at the 1 H resonant frequencies indicated in the legend, after [127]. The measurements at ν0 = 85 MHz consistently underestimate the true porosity due to signal loss from diffusion in the internal gradients.
negligible. When significant dephasing occurs as spins diffuse over a length ℓg ∼ (D0 /γ geff )1/3 , then the T2,eff measurement becomes insensitive to pore size and the relaxation time distribution is governed by details of the gradient distribution instead. Depending on the relative scaling of the lengths ℓe , ℓs , and ℓg , the observed diffusion behavior tends towards one of three limiting regimes: 1. The motional averaging (MAV) regime [123] is relevant when ℓs ≪ ℓg , ℓe . At the asymptotic limit, the magnetization decay will be governed by an exponent that varies with ℓ4s and te . This diffusive attenuation is not expected to be observed experimentally in well-connected pore networks. 2. The short time (ST) regime [124] is relevant when ℓe ≪ ℓg , ℓs ; this is equivalent to free diffusion. The magnetization decay will be governed by a surface relaxation (T2 ) exponent that varies with te and a diffusive exponent that varies with te3 . This is the optimum regime for T2 measurements as the contributions of diffusive and surface relaxation on T2,eff can be deconvolved [125]. 3. The localization (LOC) regime [126] is observed in the presence of large internal gradients where ℓg ≪ ℓe , ℓs . The asymptotic limit, where the magnetization decay is governed by a diffusion exponent that varies with te , is not expected to be observed experimentally in 2D or 3D pore geometries due to the rapid loss of signal. Outside of the asymptotic limits of these diffusion regimes, the magnetization decay is governed by a diffusion exponent that β varies with te where 1 < β < 3 and a surface relaxation exponent term that varies with te . As such, it is still possible, under certain conditions, to separate the surface relaxation and diffusion components in this pre-asymptotic regime. A thorough description of the diffusion regimes associated with internal gradients has been published in [10]. Attempts have been made to utilize diffusion in the internal gradients to determine a pore size distribution, such as the decay due to diffusion in the internal field (DDIF) experiment [128–130]. However, to be sensitive to the pore size this experiment requires a detectable degree of diffusive attenuation over the echo time te and the analysis assumes that diffusion occurs in a constant gradient. If the gradient is too small, the experiment will be dominated by surface relaxation or bulk relaxation (depending on ℓs ); if the gradient is spatially variant, the DDIF experiment will be sensitive only to features of the gradient distribution and not the pore size. As demonstrated by Mitchell et al. [10], even for a brine saturated Bentheimer sandstone with a modest magnetic susceptibility contrast measured at an intermediate field strength of B0 = 2 T, the internal gradients appear non-uniform for any practical te . Therefore, care should be taken in the application of the DDIF experiment in core analysis. The most significant consequence of internal gradients is the loss of quantification in MRI data. Several authors have considered this problem when imaging soils [131,132] and sandstones [127,133] where paramagnetic compounds are responsible for large internal gradients. Internal gradients are also a serious concern when imaging drilling mud invasion as the mud density is often increased by the inclusion of paramagnetic metal compounds [16,17]. Quantification is achieved for sandstones by measuring at very low field strengths [127,134]; see Fig. 7. For shale-free carbonates, internal gradients tend to be weak, so these can be imaged more reliably at intermediate field strengths [133,135]. Quantitative imaging methods suitable for use with carbonates and some sandstones are discussed in Section 4.2. 2.4. Multi-dimensional correlations Multi-dimensional (nD) correlations of relaxation time and/or diffusion are popular in core analysis as they enable multiple parameters to be interrogated simultaneously without a significant increase in experimental duration and provide more information than their 1D counterparts. These multi-dimensional correlations are also employed in well-logging to determine D –T2 or T1 –D –T2 correlations. Chemical or diffusive exchange can be probed using T2 –T2 or D –D measurements. Later, in Section 5.5, we discuss how these multi-dimensional correlations are being combined with profiling techniques
J. Mitchell et al. / Physics Reports (
)
–
13
in the laboratory. The mathematics behind these nD numerical inversions are challenging. Therefore, these data processing methods should only be implemented in the light of a full understanding of the limitations of the mathematics involved, which we now proceed to summarize. 2.4.1. nD inversions Correlations of relaxation and diffusion are determined by solving an nD Fredholm integral equation of the 1st kind. The method used is similar to that employed for inverting 1D data; see Section 2.1.3. The nD data are stacked into a 1D vector and compressed, allowing the inverse problem to be solved efficiently [136]. We briefly describe this process using the example of 2D NMR data. Further details are found in [92]. The NMR data are described by extending the 1D Fredholm integral equation in Eq. (6) to 2D as b (t1 , t2 ) b (0, 0)
∞
∞
K0 (u1 , u2 , t1 , t2 ) f (u1 , u2 ) d (log10 u1 ) d (log10 u2 ) + ε (t1 , t2 ) ,
= −∞
(14)
−∞
where the subscripts {1, 2} indicate the first (direct) and second (indirect) dimensions of the experiment, respectively. The data are acquired as a 2D matrix B where the direct dimension often comprises a CPMG echo train as this can be acquired in a single-shot. If a large number of echoes are acquired, this dimension can be compressed initially using a window average where the size of bins varies logarithmically [92,137]. Such window averaging requires the introduction of a precision matrix W0 to weight the compressed data appropriately. The 2D data are stacked into a 1D vector b, allowing us to solve for the ˆ The optimization problem is now stacked solution vector ˆf, which is then unwrapped to form the 2D distribution, F.
2 1/2 ˆf = arg min W0 (b − K0 f) , f≥0
(15)
2
after Eq. (7). The kernel matrix describes the expected form of the normalized data b for all discrete values in the vectors t1,2 and u1,2 1/2
1/2
as determined by the selection of experimental and processing parameters. K0 satisfies W0 b = W0 K0 f + e, where e is a stacked vector describing the experimental noise. If the two dimensions of the experiment are separable (i.e., do not share a common timebase), then the kernel function may be expressed as K0 (u1 , u2 , t1 , t2 ) = K1 (u1 , t1 ) K2 (u2 , t2 ). Therefore, the computation is made more efficient by taking advantage of the Krönecker product structure of the kernel matrix, such that K0 = K2 ⊗ K1 [136]. Now the 2D Fredholm integral equation (14) can be expressed in matrix notation as 1/2
1/2
W1 B = W1 K1 FKT2 + E,
(16)
where T denotes the transpose of a matrix, and assuming the window averaging has been applied only in the first dimension. The optimization problem becomes
1/2
Fˆ = arg min W1 F≥0
2
B − K1 FKT2 ,
(17)
2
where ∥. . .∥2 is the Frobenius norm of the matrix. To reduce the size of the optimization problem, the kernel matrix and data are compressed by obtaining the singular value decomposition of K0 and truncating the singular values to the noise floor. In this way, only the significant singular values are fitted. To improve efficiency, the singular values are determined for K1 and K2 separately (where appropriate) and used ˜ 0 ; the corresponding compressed data vector is b˜ so the reduced optimization to generate the compressed kernel matrix K problem becomes
2 ˆf = arg min b˜ − K˜ 0 f . f≥0
(18)
2
˜ 0 may have fewer rows than columns. This This reduced optimization problem is similar in form to Eq. (7), although K under-determined problem can be resolved with Tikhonov regularization. To implement Tikhonov regularization, Eq. (18) is extended with a weighting parameter λ that controls the degree of smoothing: 2 ˆf = arg min b˜ − K˜ 0 f + λ ∥f∥22 . f≥0
(19)
2
This large-scale quadratic programming problem (QPP) to Eq. (8) with L = I (identity matrix) and is solved is equivalent following the method in [138]. A fitting vector c =
˜ 0 f /λ is substituted in Eq. (19) and determined by an unconb˜ − K
strained minimization [136,138]. The optimum solution c is then interpolated to give the required distribution ˆf. Optimizing the choice of λ can be achieved using a variety of search algorithms, such as the L-curve [139], Bulter–Reeds–Dawson (BRD) [138], or GCV [93] methods. In a recent comparison of these optimization approaches, GCV was found to be the most efficient [92]. Alternative methods to Tikhonov regularization have been applied to solving the optimization problem; in well-logging a maximum-entropy-like approach is favored [137] as this method provides a rapid solution with significant smoothing, appropriate for inverting low SNR data.
14
J. Mitchell et al. / Physics Reports (
)
–
Fig. 8. Pulse sequence schematic for T1 –T2 relaxation time measurements via (a) inversion recovery [141] and (b) saturation recovery [145,146]. In all cases, the thin and thick vertical bars represent 90° and 180° rf pulses, respectively. The saturation comb of 90° pulses is repeated nc times. Homospoil gradients may be included between the pulses in the T1 encoding portion to eliminate unwanted coherence pathways. The initial magnetization at time t = 0 is set to (a) b(0) = −beq or (b) b(0) = 0. The spin ensemble is allowed to recover on the z-axis for a time τ1 before being interrogated by a CPMG echo train. The experiment is repeated for a range of τ1 times to encode T1 , followed in each case by n echoes separated in time by te = 2τ2 to encode T2 .
2.4.2. Relaxation time correlations T1 –T2 correlations have been applied to rocks [58,101,140] and oil/brine systems [141], where the additional information in the 2D distributions can simplify distinction of the oil and aqueous fluid-phase signals [134]. T1 –T2 correlations also provide the relaxation time ratio T1 /T2 which is related to the strength of surface interaction between the imbibed liquid and the solid pore matrix [142–144]. The ratio T1 /T2 is therefore useful as an indicator of wettability [143]. The 2D T1 –T2 data are generated by combining either inversion or saturation recovery measurements of T1 , see Section 2.1.1, and a CPMG measure of T2 , see Section 2.1.2; relevant pulse sequences are shown in Fig. 8. The kernel matrix K0 describing the T1 –T2 correlation is separable, with T2 in the first (direct) dimension and T1 in the second (indirect) dimension. The appropriate kernel functions are therefore K1 = exp {−2nτ2 /T2 } and K2 = 1 − c exp {−τ1 /T1 } where c will have integer values of 1 (saturation recovery) or 2 (inversion recovery). Spatially resolved T1 –T2 correlations are described in Section 5.5. Exchange times can be determined using T2 –T2 experiments [142,147], wherein a pair of CPMG echo trains are separated by a z-storage interval of duration tstore . Off-diagonal peaks are observed in the T2 –T2 distribution when spins move between environments characterized by different T2 relaxation rates on the time-scale of the storage interval. The change in intensity of the off-diagonal peaks is related to the amount of exchange. Models for diffusive exchange have been proposed, from which exchange times are calculated [147,148]. Exchange has been observed between micropores and macropores in carbonates [149] and between regions of different gradient strength in sandstones [10,122]. Exchange can also be estimated from D –D experiments [150]. At present, spatially resolved exchange measurements have not been used for core analysis. 2.4.3. Diffusion-relaxation correlations In NMR well-logging tools, the inherent magnetic field profile of the unilateral device provides a fixed field gradient (FFG) ∇ B0 that can be used to measure D . As the gradient amplitude is constant, diffusion is encoded by incrementing the echo time tse of the initial spin echo. However, this has the consequence of simultaneously encoding D and T2,eff in the signal decay. By acquiring a CPMG echo train after the diffusion encoding, the effects of relaxation are decoupled in the D –T2,eff correlation [57,101,102]. In practice, a double spin echo is acquired in the FFG-2SE experiment to provide motion compensation; the pulse sequence is shown in Fig. 9(a). As both dimensions of the experiment depend on tse , the kernel function is not separable and has the form
1 2tse + nte 3 K0 (D , T2 , tse , nte ) = exp − γ 2 g 2 D tse exp − 6 T2,eff
,
(20)
where n is the number of echoes in the CPMG train and te is the echo time separation. It is desirable to conduct complementary experiments in the laboratory on core-plugs to support results obtained downhole and improve the interpretation of the logging-tool measurements [151]. A unilateral magnetic, or a bench-top magnet with a quasi-static field gradient, could be used in the laboratory to emulate the logging-tool performance. However, this approach is limited by: 1. The inherent slice-selection of the rf pulse in the presence of the field gradient that results in a SNR reduction due to the small excitation volume. 2. Significant 19 F signal contamination will occur when fluorinated components are used in the core-holder or magnet assembly. It is preferable, therefore, to study core-plugs on the laboratory-scale using a standard magnet design that does not include a grossly inhomogeneous field profile. The pulsed field gradient double spin echo (PFG-2SE) sequence was devised to provide a close approximation to the logging-tool experiment [103], see Fig. 9(b). The diffusion part of the sequence is emulated by matching the gradient area experienced by the spins in both experiments, i.e., gtδ (rectangular gradient pulses). The kernel
J. Mitchell et al. / Physics Reports (
)
–
15
Fig. 9. Pulse sequence schematics for the equivalent diffusion-relaxation experiments (a) FFG-2SE as implemented on well-logging tools and (b) PFG-2SE as implemented in the laboratory. In the FFG-2SE sequence a background gradient of amplitude g is present at all times; the initial echo interval tse ≡ t∆ (tse /2 ≡ tδ ) is incremented in successive scans to encode diffusion. The spin echoes are narrowed by the presence of the gradient and the CPMG echo train measurements T2,eff . In the PFG-2SE sequence, pulsed gradients of amplitude g are applied to match the gradient area (gtδ ) experienced by spins in the FFG-2SE sequence. Again, tse ≡ t∆ is incremented in successive scans, along with tδ < tse . The echoes are not acquired in the presence of a gradient in the PFG-2SE sequence, so the CPMG echo train measures T2 .
Fig. 10. A schematic of the APGSTE-CPMG pulse sequence, based on the ‘‘13-interval’’ sequence of Cotts et al. combined with the standard CPMG echo train. Thin and thick vertical lines indicate 90° and 180° rf pulses, respectively. Timings on the diffusion encoding portion are specified following the notation of Tanner [95]: the echo time is tse = 2τse and the observation time (defined as the time between leading gradient edges) is t∆ . Rectangular gradient pulses are shown with duration tδ and amplitude g; the gradient amplitude is incremented in successive experiments. Homospoil gradients are indicated by gray rectangles. Data are acquired only in the CPMG portion with each train comprising n echoes separated in time by te . An extra z-storage (stimulated echo) delay of tstore is included prior to the CPMG sequence.
function for the PFG-2SE experiment is
2tse + nte tδ 2 2 2 exp − . K0 (D , T2 , tδ , tse , nte ) = exp −2γ g D tδ tse − 3
T2
(21)
The kernel is simplified if tδ can be replaced with (tse − tc ) /2 where tc is a constant time difference between tδ and tse ; this is merely a question of careful pulse sequence parameter selection. The only significant difference between the FFG-2SE and PFG-2SE experiments is found in the CPMG acquisition. In the FFG case, an effective relaxation time T2,eff is measured. However, in the PFG case, the CPMG echo train is acquired after the gradient pulses have been switched off, so – with the proviso that internal gradients are negligible – the T2 relaxation time is observed instead. Methods have been proposed to make the kernel function separable by virtue of variable substitution, although these methods require undesirable truncation or extrapolation of the data matrix [152]; it is preferable to construct the full kernel matrix when possible [103,125]. At high magnet field strengths, the APGSTE diffusion sequence is preferable as it provides some compensation against internal gradients. The 13-interval Cotts sequence is combined with the standard CPMG experiment [153] as shown in Fig. 10 [118]. An additional z-storage delay of duration tstore is incorporated between the diffusion encoding and acquisition sections of the experiment [154]. This extra delay provides time for eddy current stabilization prior to the CPMG sequence (which is sensitive to off-resonant spin dynamics [73]) and ensures the echo time in the CPMG sequence is constant for all echoes [118]. The quality of the diffusion data may also be improved by including dummy gradient pulses for eddy current stabilization at the start of the sequence [155]. Processing of the data acquired using the APGSTE-CPMG experiment (and similar pulse sequences [152]) is simplified by virtue of the two dimensions of the kernel function being separable. For the sequence shown in Fig. 10 (assuming rectangular gradient pulses and constant gradient pulse timings), the kernel functions are K1 (D , g ) =
1 2
tδ 2 2 , exp −tδ g D t∆ −
(22)
3
and K2 (T2 , nte ) =
1 2
exp −
t∆ + tstore − τse T1
−
2τse T2
exp −
nte T2
.
(23)
As the timings are fixed and known, the first exponent in K2 is included readily in the optimization problem. If the ratio of T1 /T2 is known and approximately constant over all T2 , this can be substituted to simplify Eq. (23). Generally,
16
J. Mitchell et al. / Physics Reports (
)
–
Fig. 11. Typical orientation of a core-plug in (a) a superconducting magnet and (b) a permanent magnet. The magnetic field B0 is aligned along the z-axis always; the rf B1 field is applied in the x–y plane. Common assignments of the read, phase, and slice directions are shown.
2τ (1/T2 − 1/T1 ) vanishes for bulk liquids and is ignored in rocks, and the remaining (t∆ + tstore + τ ) /T1 is assumed to be constant and therefore removed from the optimization problem by normalization of the data [152]. The factor 21 in each kernel function is associated with a stimulated echo that returns only half the initial magnetization [70]. A rf phase cycle appropriate for use with the APGSTE-CPMG sequence is given in [118]. 3. MRI in petrophysics: fundamentals The basics of MRI are covered in this section. To begin, we address some of the ambiguities that have arisen in terminology and definitions at the interface between MRI and core analysis. The theory of frequency, phase, and slice encoding methods is given in Section 3.2. 3.1. Definitions and terminology As NMR is such a diverse subject — spanning the disciplines of physics, chemistry, and medicine — it is perhaps not surprising that ambiguities have crept into the terminology and nomenclature pertaining to the wide variety of techniques available to the experimentalist. Similarly, vagaries have occurred at the interface of NMR/MRI with core analysis. We wish to be precise in the definitions and naming conventions used throughout this review to provide clarity for the reader. Taking a holistic approach, we attempt to use definitions that will be meaningful to all NMR scientists, not just those working in core analysis or well-logging. Here, we note terms that have evolved multiple definitions. All terms are defined on first use, though lists of nomenclature and acronyms are included in Appendices A and B for reference. At the most fundamental level, we use NMR to mean any magnetic resonance technique that produces data other than an image or profile, e.g., relaxation time, diffusion coefficient, chemical spectrum. We take MRI to refer to any measurement that includes a spatial dimension. Alternatives can be found in the literature, such as NMRI, although we retain the three-letter acronyms here. Within the remit of MRI, an ‘‘image’’ is assumed to have two or three dimensions, whereas a ‘‘profile’’ (or ‘‘image profile’’) has only a single dimension. In the core analysis literature, it is common to find distributions of relaxation time or diffusion coefficient referred to as ‘‘spectra’’. To avoid confusion between these data and chemical shift δ , we take the convention that a NMR spectrum is a distribution over frequencies only. We do not, after all, refer to an image profile as a spectrum even though it is a distribution over position. Similarly, the numerical inversion method used to generate relaxation time or diffusion coefficient distributions is often referred to incorrectly as an inverse Laplace transform (ILT). An inversion by the real-valued Laplace transform is a special class of first-kind Fredholm integral equation that does not apply directly to NMR data; specifically, the real valued ILT does not yield a distribution constrained to be non-negative. Accordingly, we abbreviate the concept of finding a numerical solution to a Fredholm integral equation of the first kind as a ‘‘numerical inversion’’. Chemical shift imaging (CSI) is another term that is used inconsistently. In some MRI core analysis literature, any technique that provides sensitivity to different chemical species (e.g., relaxation time or diffusion coefficient) is referred to as a CSI experiment. We prefer to reserve CSI for n-dimensional experiments where the nth dimension is a chemical spectrum or subset thereof. Alternative methods of identifying chemical species without spectroscopic resolution are referred to here as ‘‘chemically sensitive’’ MRI techniques. Another definition that requires clarification is the spatial orientation of the sample within the magnet. It is typical to define the cartesian axes x, y, and z relative to the magnetic field, i.e., B0 is always oriented along the z-axis. This convention can cause confusion as the relative orientation of B0 to a cylindrical core-plug will change depending on magnet geometry. For example, a superconducting magnet has B0 aligned along the magnet bore, so z will typically be parallel to the long axis of a cylindrical core-plug. A low-field permanent magnet will typically have B0 aligned perpendicular to the magnet bore, so z will cut across the long axis of cylindrical core-plug, with y (or x) parallel to the long axis. The rf field B1 is always applied in the x–y plane, perpendicular to B0 . The main orientations are summarized in Fig. 11. We now define the concepts of porosity and saturation for those unfamiliar with the scaling laws and units often employed for these values. Porosity φ of a rock is defined either as a fraction of the total plug volume V , or in porosity units (p.u.) on a scale of 0 to 100 where 1 p.u. is equivalent to φ = 0.01 porosity. Obviously, porosity can be stated as an
J. Mitchell et al. / Physics Reports (
)
–
17
absolute volume in cm3 , equivalent to the dimension ‘‘pore volume’’ (PV) such that 1 PV ≡ φ V . Similarly, fluid saturation S is quoted as either a fraction of the total porosity, or in saturation units (s.u.) on a scale of 0 to 100, where 100 s.u. means the entire pore volume is filled with that fluid. Saturation of a particular fluid-phase is denoted by an appropriate subscript, such as So for oil saturation or Sw for water (brine or aqueous fluid-phase) saturation; remaining or residual oil saturation is indicated by an additional subscript ‘‘r’’, e.g., Sor . For an oil recovery experiment, the state of the flood is indicated by a (initial) superscript term; for example, the initial oil saturation would be quoted as So , whereas the final oil saturation might (final) be given as Sor . Finally, we note that in the literature it is common to see both porosity and saturation quoted in %; this scaling can be confusing without a clearly defined reference volume. 3.2. Basic imaging theory Spatial position is encoded in the reciprocal k-space signal b(k) and the profile amplitude bˆ (r) is obtained from the FT of b(k); k-space is therefore the Fourier inverse of real space [5]. The Larmor frequency ω of the spins at position r is related to the applied magnetic field gradient vector g by (in the laboratory reference frame)
ω(r) = γ B0 + γ g · r.
(24)
When k-space data are acquired on resonance in the rotating frame, the voxel intensity is obtained by the 3D FT of b(k) such that
bˆ (r) =
b(k) exp [+i2π k · r] dk.
(25)
3.2.1. Frequency encoding In medical MRI it is typical to acquire a 1D profile (or single dimension in multi-dimensional imaging) using a frequency encoded gradient echo [156]. As the amplitude of a gradient echo depends on T2∗ , it is desirable when imaging porous media to include a coincident frequency encoded spin echo (dependent on T2 ); this modification typically requires only the addition of a refocusing rf pulse [156]. Frequency encoding has the advantage of providing a profile in a single scan where the resonant frequency of the spins is encoded spatially via the application of a constant amplitude read gradient. A frequency encoded spin echo provides a single line of k-space points. Frequency encoded profiles are not used for quantitative measurements in intermediate-field imaging because the intensity of the profile is weighted by both T2 and T2∗ relaxation. A particular complication of frequency encoding is the convolution of the spectral linewidth (Fourier inverse of T2∗ ) and the profile. However, quantitative frequency encoded profiles have been obtained in very low-field experiments [134]. A frequency encoded spin-echo profile is obtained by applying a magnetic field ‘‘read’’ gradient such that excited nuclear spins precess with a frequency determined by their position. If the magnetic field gradient is applied for a time tr , then the range of wavenumbers spanned is given by
± k(max) = ±
γ gtr , 4π
(26)
where kmax is the maximum deviation from k = 0 m−1 , the central k-space position corresponding to ω0 . Each point in k-space will have a width defined by
∆k =
γ g ∆t , 2π
(27)
where ∆t is the dwell time (the time over which the continuous analogue rf signal is summed to form a discrete digital datum). The field of view (FOV) of the profile is then determined by the reciprocal relation FOV = 1/∆k such that the spatial resolution is ∆z = FOV/m where m is the number of complex data acquired. Modern digital spectrometers have the capability to over-sample in the time domain, so the practical FOV will be defined by the filter bandwidth fBW rather than the reciprocal dwell time. Over-sampling permits the removal of the noise contribution outside the filter bandwidth. Frequency encoded profiles must be acquired with sufficient bandwidth to exceed the frequency range
∆ν0 =
γ gL , 2π
(28)
where L is the sample length. This increased bandwidth reduces the SNR per scan compared to a bulk sample measurement (where a narrow band filter may be employed) by allowing more noise into the receiver. The frequency bandwidth is of particular significance in low-field measurements but less so in intermediate to high-field experiments where the SNR is higher anyway. Typically, a frequency encoded profile can be measured only in a single direction (either x, y, or z, or linear combination thereof) in any one experiment. Read profiles acquired in multiple directions can be reconstructed using back-projection algorithms [157]. In modern MRI, it is far more common to use FT imaging reconstruction by combining frequency, phase, and slice techniques [5,156]. Complicated k-space trajectories may be defined using time-varying read
18
J. Mitchell et al. / Physics Reports (
)
–
gradients, although such acquisition techniques have not yet been applied in core analysis. Rather, a single read direction is usually chosen to correspond to the long axis of the core-plug, see Fig. 11, as this provides the most efficient method for acquiring the largest number of k-space points (pixels) to achieve isotropic image resolution. The total intensity of the frequency encoded read profile is determined by the signal acquired at k = 0 m−1 (the center of k-space). However, due to the combination of g and ∆t in Eq. (27), it can be difficult to acquire a k-space datum that corresponds to, or at least averages uniformly about, k = 0 m−1 . Improved profile quality can be achieved by careful experiment design: ensuring a datum is acquired that provides a good estimate of the signal intensity at the center of k-space. The use of read gradients is generally impractical when samples exhibit very short T2∗ relaxation times as the rapid FID decay masks any spatial information in the signal. However, imaging of such broad line samples can be achieved using very strong read gradients that dominate the inherent linewidth of the material. Suitably strong (static) gradients are found in the stray field of superconducting magnets. By positioning the rf resonator and sample in a region of the stray field where the gradient is both large and constant over the volume of interest, a 1D profile is obtained. This technique is known as stray field imaging (STRAFI) [158]. STRAFI has been used to monitor water ingress into sandstone plugs [159], although its application to rocks has been limited; several drawbacks of the STRAFI technique mean it is not appropriate for core analysis. A STRAFI profile can only be acquired in 1D (due to the alignment of the static magnetic field gradient) and the sample must be rotated to obtain 2 or 3D images. Furthermore, the sample size is limited to the uniform region of magnetic field gradient (typically on the order of millimeters) and 19 F used in core-holders will be excited following the arguments rehearsed in Section 2.4.3. Broad line imaging techniques using phase encoded pulse gradients are more appropriate for core analysis, as discussed in Section 4.2. 3.2.2. Phase encoding To select a k-space position using phase encoding, a ‘‘phase’’ gradient of variable amplitude g is applied for a time tp . T2∗ blurring is negated when using phase encoding because the acquisition time (i.e., time after initial excitation) of each datum is fixed. The k-space point sampled at a particular gradient strength is determined by k=
γ gtp . 2π
(29)
By adjusting the gradient amplitude between ±g (max) over m subsequent scans, a line of k-space is traversed. The corresponding FOV is FOV =
mπ
γ g (max) tp
,
(30)
and the image resolution is
∆y =
π γ g (max) t
.
(31)
p
The phase encode gradient can be ramped in multiple directions (x, y, or z) in subsequent experiments, so that k-space may be traversed in 1, 2, or 3D. When all the image dimensions are acquired in this manner, the technique is referred to as a pure phase encode method. 3.2.3. Slice selection A spatial slice is selected by applying a shaped rf pulse in the presence of a gradient. Slice selection differs from frequency and phase encoding in that the obtained signal intensity is directly proportional to the spin density in the slice. A profile may therefore be obtained by moving the location x0 of the resonant slice sequentially through the sample; no further data processing is required. The width of the slice ∆xs is determined by the frequency bandwidth ∆ν of the rf pulse and the gradient amplitude g as
∆ xs =
2π ∆ν
γg
,
(32)
where the rf pulse bandwidth is a property of the pulse shape. The position of the slice is determined by the frequency offset, γ gx0 /2π . The inverse relationship between frequency and time means a narrow slice width ∆xs is achieved only with a long rf pulse duration, often ≫ 100 µs when using a modest gradient amplitude. The slice shape is the Fourier inverse of the rf pulse shape, so an ideal boxcar slice is obtained from a sinc pulse with an infinite number of oscillations. The practical quality of the slice will be influenced by the truncation of the rf pulse width to a realistic duration, coupled with the ability of the spectrometer and rf amplifier to reproduce the required shape. It is often easier to define a Gaussian rf pulse and accept that the slice shape will also be a Gaussian; the simplified waveform provides better reproducibility and control for the experimentalist. The slice selection gradient of amplitude g is applied for a time ts , ideally equal to the duration of the rf pulse. During this time, the excited portion of the spin ensemble will acquire a phase shift. It is therefore necessary to apply a rephasing gradient lobe of approximate area gts /2 after the slice select gradient to remove the phase shift. The exact
J. Mitchell et al. / Physics Reports (
)
–
19
area of the rephase gradient lobe relative to the slice select gradient is dependent on details of the rf pulse shape [156]. Slice selection can be inappropriate when imaging short T2∗ samples due to the long rf pulse durations. Physical slice selection can also be achieved by shielding parts of the sample from the rf field with metal plates. Lucas et al. used this approach to image sections of a long rock core [160]. This approach, although impractical for widespread application, is most useful for preventing signal detection from regions of sample outside the rf coil, or masking entry effects in flow experiments. However, most modern spectrometers possess the capability to achieve the same selective excitation with a shaped rf pulse in the presence of a weak slice gradient [161]. This latter method is preferable as the experimentalist has more control over the quality of the rf pulse: the presence of a metal shield will distort the rf field. Notwithstanding, physical slice selection has the potential to be useful when running pulse sequences that contain a large number of closely spaced rf pulses, such as CPMG, where the implementation of a soft rf pulse and slice gradient on each excitation is not practical. An alternative to a shaped rf (soft) pulse is DANTE-Z slice excitation [162–165] which employs a series of hard rf pulses to saturate unwanted (out-of-slice) spins and leave the required (in-slice) spins on the z-axis. DANTEZ slice selection can be followed by any imaging sequence, including CPMG-like echo trains as it is insensitive to spin phase, unlike conventional soft rf pulses that can introduce phase artifacts when preceding multi-echo imaging techniques [166]. However, this technique works only when T1 ≫ T2 as out-of-slice spins saturated initially will recover and be detected following imperfect 180° rf pulses. 4. High resolution MRI in core analysis The primary application of MRI has been in medicine, where anatomical images of humans and rodents are used for clinical diagnosis and biomedical studies, respectively. The concept of applying MRI to medical diagnosis was suggested by Damadian [167]. The practical implementation of MRI was demonstrated by Mansfield and Grannell [168], and Lauterbur [169]. The first full-body MRI scanners were constructed nearly a decade later [51]. Since then, there has been a drive in magnet technology to attain higher field strengths for improved SNR (leading to improved image resolution) and spectroscopic resolution for fMRI [156]. High resolution images of fluid-saturated rocks are obtained at intermediate to high field strengths as will be discussed in Section 4.1. Much of the work in the literature has focused on mass transport, such as forced displacement of oil by water, or drainage against capillary forces. However, the lack of quantification when using fast medical MRI protocols has meant that these methods have lost favor in the core analysis community. Other imaging techniques are perceived to perform better for core analysis applications [12]. These quantitative methods are discussed in Section 4.2. Notwithstanding, intermediatefield images are an excellent non-destructive and non-invasive measure of sub-millimeter-scale heterogeneity, allowing the identification of structural irregularities in short core-plugs. In this section we review the main imaging methods and consider how they have been applied in core analysis. All gradient orientations are given for core-plug imaging in superconducting magnets according to Fig. 11(a). 4.1. Heterogeneity mapping — rapid imaging techniques Imaging of materials (soft solids, porous media, solids, etc.) has proven more difficult than biomedical systems due to the short relaxation times associated with these samples. The earliest example of MRI in core analysis was presented by Gummerson et al. who profiled 1D spontaneous ingress of water into sandstones, a limestone, and cement based materials using a low-field magnet [15]. Half a decade later, 2D images of rock core-plugs were published by Rothwell and Vinegar [27] with a description of their methodology given in [170]. A comparison of MRI and X-ray CT was presented in [171]. Fast imaging techniques have also been used to monitor gel formation in rocks [190]. Early 2D imaging was achieved with projection reconstruction [27,172,173]; however, this modality is implemented rarely nowadays and so we have chosen not to discuss the details here. Rather, we focus on more common and versatile Fourier imaging methods. High resolution 3D imaging is possible with the spin-warp pulse sequence [174,175]. The standard implementation utilizes a frequency encoded echo in the read dimension, combined with two phase encode gradients to achieve full 3D k-space sampling. There is inherent relaxation contrast in the images generated by the spin warp method, although the timing of each k-space acquisition is constant, resulting in the final image having a uniform weighting within the limitations of the frequency encoded spin echo acquisition (i.e., T2∗ ). The spin warp sequence does not provide the temporal resolution required for medical imaging, although it is suitable for heterogeneity mapping in static core-plugs. The basic spin warp pulse sequence is shown in Fig. 12. A full isotropic image of a core-plug will take several hours using the spin warp sequence. For example, acquiring a spin warp image of typical size 128 × 64 × 64 voxels (read × phase × phase) where each scan has a duration texp = 10 s (including recovery delay), and 4 repeat scans are averaged, will take in excess of 45 h. The inherent relaxation time weighting in spin warp images has been used to identify voxels of differing mean porosity [176]. An example 3D spin warp image of a core-plug is shown in Fig. 13. Although many authors have applied spin warp imaging to determine rock heterogeneity [8,177–185], faster imaging protocols are preferable as they can be extended more readily to time-varying systems. Fast imaging techniques are preferred for clinical studies to minimize the time the patient spends in the magnet. Such protocols are useful in core analysis for investigating time dependent systems, as explored in Section 4.3. A popular technique
20
J. Mitchell et al. / Physics Reports (
)
–
Fig. 12. 3D spin warp imaging sequence utilizing one read gradient and two phase encode dimensions. The thin and thick vertical lines represent 90° and 180° rf pulses, respectively. The spin echo is centered at time te after the initial rf excitation pulse. The phase encode gradients are incremented individually in successive experiments.
Fig. 13. 3D reconstruction of a spin warp image of a brine saturated Ninian sandstone core-plug. The full image contained 256 × 128 × 128 voxels with an isotropic pixel resolution of 0.3 mm in each direction. A portion of the image is cut away to reveal the shale band (dark line) traversing the plug. The gray scale ranges from black (no signal) to white (maximum signal). Source: Image reproduced with permission from The Royal Society [180].
in clinical imaging is echo planar imaging (EPI) [186,187] and its various derivatives. However, as EPI is a gradient echo imaging sequence, it is not well suited to rocks with short T2∗ relaxation times. Furthermore, EPI images are susceptible to artifacts from inhomogeneities in B0 and B1 , and chemical shift effects [188]. The hybrid π -echo planar imaging (PEPI) — originally called ‘‘180° EPI’’ — reduces the sensitivity to off-resonance effects but relies on good quality 180° rf pulses. PEPI was developed for studies of materials [24,189], relying on spin echoes (and hence T2 relaxation) rather than gradient echoes. A modern version of the PEPI pulse sequence is shown in Fig. 14. Like the spin warp pulse sequence in Fig. 12, PEPI utilizes frequency encoded spin echoes in the read dimension and phase encode gradients in the other two dimensions. However, PEPI offers a temporal improvement over the spin warp image acquisition because a 2D k-space raster is acquired on each scan. For example, acquiring a PEPI image of size 128 × 64 × 64 voxels (read × phase × phase) where each scan has a duration texp = 10 s (including recovery delay), 4 repeat scans are averaged, and assuming 64 echoes are acquired (i.e., full 2D k-space sampling per scan), will take less than 1 h. Through the mid 1990s, PEPI was demonstrated as a suitable imaging method for heterogeneity mapping of core plugs [24,25]; however, it does have limitations. Examples of early images are shown in Fig. 15. PEPI offers a significant time saving over spin warp imaging but contains more T2 relaxation weighting due to the use of a CPMG-like echo train. The T2 weighting is distributed non-uniformly through the image and depends on the acquisition order of the k-space raster. Furthermore, PEPI is notoriously difficult to implement well, with imaging artifacts arising from non-ideal rf pulses and mismatched gradients. An alternative is the rapid acquisition with relaxation enhancement (RARE) sequence [26] shown in Fig. 16. Like PEPI, RARE makes use of a CPMG-like echo train to acquire multiple k-space lines in a single scan. The number of echoes n acquired is called the RARE-factor. Ideally, the RARE-factor is sufficient to encode an entire 2D image in a single scan, although this will depend on T2 . RARE is typically implemented as a multi-slice sequence for pseudo-3D imaging, which means that signal is detected from one slice while the previously excited slice recovers. In this way, the requirement that the delay between scans be 5 × T1 is circumvented. However, there are also limitations with the RARE sequence. The speed of a RARE acquisition is limited typically by the gradient (hardware) duty cycle. Soft rf pulses and slice gradients are used for the excitation and every refocusing rf pulse in the echo train. The use of soft rf pulses does limit the minimum echo time, although the slice
J. Mitchell et al. / Physics Reports (
)
–
21
Fig. 14. The PEPI-hybrid imaging sequence allows the acquisition of fast images in short T2∗ materials, such as rocks. Signal acquisition begins after the first echo; the first (unobserved) spin echo is not shown. The thin and thick vertical lines represent 90° and 180° rf pulses, respectively. The acquisition sequence is repeated n/2 times with a blipped phase gradient applied on the first phase (y) axis to step through a 2D k-space raster during the echo train. The third k-space dimension is sampled by ramping a phase gradient on the second phase (x) axis.
Fig. 15. 2D image slices extracted from 3D PEPI images of (a) a vuggy dolostone plug, (b) a partially saturated Berea sandstone plug, and (c) a Brent oilfield plug (as cored, stored under brine). Each 3D image contained 128 × 64 × 64 voxels with an isotropic pixel resolution of 1.5 mm. The gray scale ranges from black (no signal) to white (maximum signal). Images reproduced with permission from Elsevier [24].
Fig. 16. RARE pulse sequence schematic. The low and high soft (Gaussian) rf pulses provide nominal tip angles of 90° and 180°, respectively. A series of n spin echoes are acquired, where n is referred to as the RARE-factor. The phase (y) gradient is incremented on each echo. Slice select (x) gradients are applied to coincide with each soft rf pulse. The gray rectangles are homospoil gradients. The sequence can be modified to acquire T2 maps (MSME) by incrementing the phase (y) gradient on every scan rather than every echo.
gradients also act as spoilers that prevent unwanted spin coherences. The echo train generated by the rf pulses, intentionally set to provide incomplete refocusing of the spin ensemble, contains a mixture of spin echoes and stimulated echoes which are acquired superimposed; the relaxation weighting is then a mix of T2 and T1 which has advantages: in particular, the addition of T1 weighting means the echo train can be acquired over a longer time (i.e., more echoes) than in, say, PEPI that relies only on T2 relaxation. However, the mix of T1 and T2 weighting also ensures quantitative image intensities cannot be determined with the standard RARE sequence.
22
J. Mitchell et al. / Physics Reports (
)
–
Fig. 17. 2D longitudinal (left) and transverse (right) images extracted from 3D multi-slice RARE images of four brine saturated microporous limestone core-plugs (top to bottom). On visual inspection, all four plugs appeared identical; high-resolution MRI reveals subsurface structural variations. Each image has an isostatic pixel resolution of 0.5 mm. The slice thickness was 1 mm. The intensity (gray-scale) is a qualitative indication of local porosity.
Due to the nature of the acquisition, RARE images are considered qualitative and cannot be corrected easily for relaxation weighting. Nevertheless, the speed advantage is significant. For example, acquiring a multi-slice RARE image of size 128 × 64 × 64 voxels (read × phase × slice) where each scan has a duration texp = 3 s, 4 repeat scans are averaged, and assuming a RARE-factor of 64, will take approximately 15 min. Example RARE images obtained at B0 = 2 T are shown in Fig. 17. Other combinations of image acquisition protocol are achieved by combining the concepts of frequency, phase, and slice encoding [156]. Spin warp, PEPI, and RARE remain the main sequences used for rapid, high-resolution imaging of materials. Geostatistical structure analysis [191] is applied to MRI images in the form of variograms. The variogram function Υ (r) describes the correlation between Nv pairs of points separated by a vector r as 2Υ (r) =
Nv 1
Nv j=1
2 Ξ ξj − Ξ ξj + |r| ,
(33)
where the parameter Ξ is determined at position ξj . Ξ can be any property of the pixels, such as intensity (scaled to porosity), or relaxation time or diffusion coefficient if appropriate maps are acquired. By calculating variograms for different orientations of r in an image, ordered structuring in the rock is detected. Bansal et al. applied variogram analysis to a selection of core-plugs, including Berea sandstones and Indiana limestone [192]. A modern approach is given in [193]. The fast imaging techniques discussed so far have all been examples of qualitative imaging protocols. These methods of qualitative imaging are most useful when combined with the suite of other techniques available to form a coherent core analysis package [194–196]. A combination of low field relaxation time measurements and imaging has been used to characterize carbonates [197,198], leading to a method of monitoring wettability changes [199]. X-ray CT has been used to support MRI data, leading to improved rock characterization [200,201]. Padhy et al. performed MRI as part of a core analysis study of vuggy carbonates, with pore size distributions estimated by DDIF [202,203]. 3D images of partially drained rocks have also been combined with impedance measurements to determine improved resistivity indexes [204]. 4.2. Heterogeneity mapping — quantitative imaging techniques Quantitative MRI is achievable if the relaxation weighting associated with each voxel is determined. Rothwell and Vinegar [27] implemented a multi-echo back-projection imaging sequence to generate a 2D T2 map. Chen et al. attempted to correct signal loss in images of multi-fluid-phase flow by modeling the variation in T2 as a function of saturation [205,206]; they
J. Mitchell et al. / Physics Reports (
)
–
23
Fig. 18. Pulse sequence schematic for multi-echo, spatially resolved T2 mapping, after [28]. The thin and thick vertical lines represent 90° and 180° rf pulses, respectively. A series of n CPMG-like echoes are acquired in the presence of a read gradient of amplitude g and duration tr . A T2 decay is associated with each k-space point (and hence pixel after Fourier transformation). The read gradient can be applied on any one axis; multiple directions must be scanned in separate experiments. Acquiring nD T2 maps requires the addition of phase encode gradients. The addition of bipolar phase encode gradients, as in RARE (see Fig. 16), would be appropriate for extending the multi-echo sequence to 2 or 3D.
Fig. 19. Pulse sequence schematics for (a) SPI and (b) SPRITE. The rf excitation in (a) is via a 90° pulse, and in (b), via a small tip angle α pulse (where α ≪ 90° such that sin(α) ≈ α ). The usual implementation involves the acquisition of a single point at time tp after the rf pulse. Multiple points may be acquired and adjusted to a constant FOV using a z-chirp transform. The phase gradients are ramped independently between ±g (max) in either 1, 2, or 3 dimensions. For (b) SPRITE, the gradients are incremented on each loop to sample n k-space positions in a single scan.
achieved limited success as their bulk measurement of T2 did not take into account local variations in relaxation time. Majors et al. demonstrated a 1D pulsed gradient version of the multi-echo imaging sequence [28], illustrated in Fig. 18. A series of spin echoes, each acquired in the presence of a read gradient, provides a CPMG decay for each pixel in the profile. These decay curves are fitted using the method described in Section 2.1.3 to obtain b(0), the projected signal at zero time, and a T2 distribution. Analogous to the read dimension of the PEPI sequence shown in Fig. 14, this 1D multi-echo sequence is often assumed to be derived from RARE; however, since soft rf pulses are used in RARE, the spin dynamics are very different. Further discussion on the 1D multi-echo sequence is given in Section 5. The multi-slice multi-echo (MSME) sequence [156] offers a more sophisticated method of obtaining relaxation time maps in 2 or 3D. This sequence is identical in design to RARE, see Fig. 16, except the phase encode gradient is ramped over subsequent scans rather than subsequent echoes. In this way, the data include an echo train that decays with a mixture of T1 and T2 associated with each voxel. Relaxation time maps are again obtained by fitting the decay curves. A pure T2 relaxation time measurement is achieved by adding additional spoiler gradients to shift the stimulated echo contribution out of the acquisition window; note that this modification reduces the acquired signal intensity. Porosity and permeability maps have been generated from such relaxation time maps [207,208]. The MSME sequence has a substantially longer total experimental duration compared to a RARE image of equivalent spatial resolution. For example, acquiring a MSME image of size 128 × 64 × 64 voxels (read × phase × slice) where each scan has a duration texp = 3 s, and 4 repeat scans are averaged, will take approximately 14 h. Therefore, the MSME sequence is implemented typically to acquire a few slices of interest rather than a full 3D image. To ensure relaxation-corrected images are truly quantitative, some consideration should be given to macroscopic spatial variation in the rf field B1 which can lead to variable spin dynamics throughout sample [209]. Relaxation mapping can be appended to any fast imaging sequence by pre-conditioning rf pulses [210,211] (such as inversion recovery for T1 or CPMG for T2 ) on the front of the imaging sequence [156]. The image acquisition is repeated for different relaxation delays, so that the voxel intensity decays with T1 or T2 over a series of images. Relaxation time distributions can again be obtained as described in Section 2.1.3. If the relaxation rate is slow relative to the image acquisition time, this approach may be used to generate quantitative images from fast imaging protocols such as RARE [212]. Another method of obtaining quantitative MRI data is single point imaging (SPI) [30]. Basic SPI techniques are robust to hardware instabilities and sample-to-sample variability, and are therefore considered reliable when compared to other imaging protocols. These properties have made SPI-based methods popular for core-analysis applications. In SPI, each k-space point is acquired at a fixed time tp after a rf excitation pulse, in the presence of a phase encode gradient; see Fig. 19(a). Accordingly, the T2∗ relaxation weighting is constant across k-space and therefore does not distort the profile. If the acquisition time tp of each k-space point is short compared to T2∗ , then the resultant profile is considered quantitative. As such, SPI is well suited to the measurement of saturation states in rock core-plugs. It has been shown that the FID decay obtained from liquid saturated rocks is mono-exponential over a wide range of magnetic fields and rock types [10,133,213]. Therefore, calculating the absolute pixel intensity in a SPI image is straight-forward. SPI is time consuming because each k-space point is acquired in a separate scan. Accordingly, SPI has been improved in the form of single point ramped imaging with T1 enhancement (SPRITE) [31] which uses repeated applications of small tip angle rf pulses to scan a
24
J. Mitchell et al. / Physics Reports (
)
–
Fig. 20. 2D image slices extracted from 3D conical SPRITE images of a partially saturated Berea sandstone plug. Longitudinal (left) and transverse (right) slices are shown for total water saturations of (top) Sw = 26.3% and (bottom) Sw = 30.9%. The 3D SPRITE images contained 64 × 64 × 64 voxels with an isotropic resolution of approximately 1.094 mm. The phase encode time was tp = 40 µs. An entire 3D image was acquired in 2.5 min. Source: Images reproduced with permission from the American Institute of Physics [56].
complete k-space raster in a single scan; see Fig. 19(b). This pulse sequence relies on short T1 relaxation times to continually restore a fraction of the spin ensemble to the z-axis, thus providing signal over a long (∼1 s) acquisition window. SPRITE offers a significant reduction in experiment duration compared to SPI, although the required high gradient duty cycle does mean commercial hardware struggles to generate the necessary gradient waveforms that can last for several seconds. Similarly, commercial hardware may fail to deliver suitably short rf probe response times (probe dead time); ideally, the time between rf excitation and detection should be tp ∼ 10 µs. As SPI and SPRITE data are acquired in the presence of a phase encode gradient, the frequency bandwidth limit given in Eq. (28) applies. SPRITE has been used successfully to image rocks across a wide range of magnetic field strengths from B0 = 0.18 T [214] to B0 = 7 T [215–217]. However, application to very low magnetic field strengths is restricted by the inherent poor SNR arising from the small tip angle rf pulses and wide frequency bandwidth. SPRITE is inherently insensitive to multiple fluid-phases due to the lack of relaxation time weighting in the images. Methods for introducing fluid-phase contrast into SPRITE through pre-conditioning have been explored [218], although these techniques are not straight-forward to implement, nor do they allow multiple fluid-phases to be observed simultaneously with a high temporal resolution. SPRITE imaging has been used in various petrophysical studies to identify structures in core-plugs [215,219,220]. The pure phase encode nature of SPRITE lends itself readily to k-space sampling schemes other than rectilinear grids. The majority of the intensity in an image is acquired at low k positions. Therefore, under-sampling of high k positions reduces both the image fidelity and acquisition time but does not reduce the quantitative nature of the image. In core analysis, the determination of accurate signal intensity and relaxation time or diffusion coefficient is often more important than a well defined image. This bias towards quantity over quality in materials science and engineering is largely the opposite of the requirements for medical MRI. SPRITE k-space data may therefore be acquired with a spiral (2D) or conical (3D) sampling scheme [221]; see, for example, Fig. 20. SPI-based sequences are well suited to the random k-space sampling required for compressed sensing reconstructions [222] (see Section 6.1). The SNR of SPRITE is improved by sampling the origin of k-space first, as encompassed in the 1D centric (double half-k) SPRITE acquisition, where the phase encode gradients are ramped from g = 0 to +g (max) in one scan, and g = 0 to −g (max) in the subsequent scan. When using complicated sampling schemes (such as 3D cones), the origin of k-space may be sampled multiple times, providing an improved SNR on the total image intensity. Furthermore, multiple points may be sampled after each rf excitation in both SPI and SPRITE, as indicated in Fig. 19. However, as each point has a different associated tp , a series
J. Mitchell et al. / Physics Reports (
)
–
25
Fig. 21. Pulse sequence schematics for (a) the original SESPI sequence and (b) the modified SESPI sequence. The thin and thick vertical bars represent 90° and 180° rf pulses, respectively. The phase gradients are ramped independently in either 1, 2, or 3 dimensions. In (a) the bipolar phase encode gradients may be incremented on each echo to sample a line of k-space in a single scan (turbo-SPI), or incremented on sequential scans to provide a T2 map. The modified SESPI sequence (b) only provides a T2 map.
of images — each with a different field of view determined by tp according to Eq. (30) — are obtained. The fields of view are standardized using the z-chirp transform [223,224]. SNR can also be improved by customizing the digital filter bandwidth used to sample each k-space location [225], although implementation of this method can be impractical depending on the hardware available. When imaging a rock with T2 ≫ T2∗ , improved SNR is achieved with the spin echo single point imaging (SESPI) sequence [32,226], shown in Fig. 21(a). Bipolar phase encoding gradients of amplitude ranging over ±g (max) , each applied for a time tp , are located on either side of every echo. This pulse sequence is equivalent to the phase encode scheme used in RARE. Only the echo amplitude is acquired in the SESPI experiment to provide a single k-space point. As the signal is not acquired in the presence of a gradient, there is no FOV limitation on the minimum frequency bandwidth that can be used. A narrow band digital filter provides an additional SNR improvement over SPI. The advantage of this SESPI sequence is the simplicity of the spin dynamics: the spin ensemble is coherent during the rf transmit and receive intervals. However, imaging artifacts will occur if the positive and negative gradient lobes are not matched precisely. Like SPI, the phase encode gradients may be applied on all axes to generate a 2 or 3D image. If the phase gradients are incremented on every echo, a line of k-space is traversed in a single scan; this acquisition mode is referred to as turbo-SPI [226,227]. Alternatively, the phase encode gradients may be incremented on each scan to provide a CPMG-like decay, and hence T2 distribution, for each voxel [32]. 2D T2 mapping by this method is discussed in [222]. A modification to the SESPI sequence offers a good compromise when a high SNR is required from a short T2 sample [228]. The modified SESPI sequence, Fig. 21(b), has only a single phase encode gradient of variable amplitude ranging over ±g (max) , applied for a time tp , between the leading excitation and refocusing rf pulses. This sequence is derived from the second phase encode dimension in PEPI. The position of the gradient pulse allows the subsequent echo amplitudes in the CPMG-like echo train to be acquired with a reduced echo time te′ = 2τ ′ . However, this sequence has two limitations: (1) the initial echo time te = 2τ is still long so very short T2 components will go unobserved, and (2) the spin dynamics are complicated. Despite the long initial echo time, the subsequent echoes, closely spaced, will provide the maximum achievable signal from short T2 components. To overcome the phase shifts that accumulate by repeated reversal of the incoherent spin ensemble, the xy-16 phase cycle [229] is applied to the 180° rf pulses along the echo train. Consequently, only the decay described by every 16th echo is a good approximation of the T2 relaxation behavior. Therefore, the effective echo time becomes te′ = 32τ ′ if accurate T2 relaxation rates are to be determined. Additionally, composite rf pulses have been used to improve the quality of the echo shape [228]. These long rf pulses further limit the minimum echo time achievable. Notwithstanding, the direct implementation, as shown in Fig. 21(b), provides a reasonable T2 decay when fitting all the echo amplitudes. We note that the initial echo time of the modified SESPI sequence is still short when compared to say, frequency encoded echoes, as the phase encode method relies only on the total gradient area, not the gradient shape. Gradient switching and stabilization delays can therefore be minimized in the modified SESPI experiment. The improved SNR of the SESPI sequences makes them suitable for implementation at very low magnetic field strengths. Further discussion on these imaging sequences in 1D, along with examples of T2 maps and applications, is given in Section 5. 4.3. Fluid transport A particular advantage of MRI for core analysis is the capability to spatially resolve the distribution of multiple fluid-phases within a rock core-plug. Chemical spectroscopy provides the best fluid-phase discrimination and allows simultaneous detection of multiple liquid components in core analysis. However, obtaining high-resolution chemical spectra is often not possible due to line broadening caused by the magnetic susceptibility contrast in the heterogeneous samples. There are two distinct methods of chemical shift imaging (CSI), defined here as (1) spatially resolved chemical spectra (i.e., high resolution spectra and a low resolution image), or (2) chemically resolved images (i.e., high resolution images and a low resolution spectrum or chemically selective detection). We now consider these two approaches to CSI in turn. Spatially resolved spectra are obtained using either pure phase encoding [230] or volume (slice) selection [231]. Both methods rely on the acquisition of a FID in the absence of an imaging gradient. Example pulse sequences are given in Fig. 22.
26
J. Mitchell et al. / Physics Reports (
)
–
Fig. 22. Pulse sequence schematics for (a) 3D phase encoded CSI and (b) 3D volume selective CSI. In (a), the thin and thick lines indicate 90° and 180° rf pulses, respectively. The gray rectangles are homospoil gradients that can be applied on any axis. In (b), the short and tall Gaussian waveforms indicate 90° and 180° soft rf pulses, respectively. The phase and slice methods may be combined as demonstrated in [232].
Phase encoding and slice selection may be combined as required; Dereppe et al. acquired 3D oil and water images from coreplugs using 1D slice selection and 2D phase encoding [232,233]. Regardless of the method for spatial selection, the FIDs are Fourier transformed to provide a chemical spectrum for each pixel. The area under the peaks of interest are integrated (or determined by peak fitting) to provide chemically selective images. This analysis technique is robust even when significant overlap of the spectral peaks occurs, as is common for liquid-saturated rocks due to the inherent line broadening caused by the inhomogeneous sample. When the spectral lines are well separated, a soft rf pulse (in the absence of a slice selection gradient) is used to selectively excite the spins in one chemical species. Dechter et al. used a chemically selective rf pulse to saturate one component in an oil/water saturated rock, thereby obtaining signal only from the other component [234,235]. The use of a chemically selective pulse is advantageous for dynamic studies as an image of a single fluid-phase is obtained directly. However, good selection is hard to achieve in liquid-saturated rocks, again due to the broad and therefore overlapping lines. Further, the soft rf pulse required to achieve good chemical selectivity will have a long duration, which will be unsuitable for rocks with short T2 or T1 relaxation times. Chemical selectivity can also be achieved by incorporating an additional phase evolution delay in the standard spin-warp imaging sequence [236]. This technique relies on the chemical sensitivity inherent in the FID to detect different species selectively. An incremental phase evolution delay tϕ is added between the excitation rf pulse and the read gradient in the spin-warp sequence [237]. If Nϕ evolution times are used, tϕ,n =
n 2 Nϕ − 1 (ν1 − ν2 )
,
n ∈ 0, . . . , Nϕ − 1 ,
(34)
where (ν1 − ν2 ) is the difference between the resonant frequencies of the chemical components (e.g., oil and water). Dixon demonstrated that two chemical species, whose central frequencies (ν1 and ν2 ) are well separated, can be distinguished with Nϕ = 2 [236]. However, in heterogeneous rock samples, phase artifacts are common due to distortions in the background magnetic field. Horsfield et al. demonstrated that background phase maps could be inferred and used to correct the chemically selective images; this processing required Nϕ = 8 [238,239]. An improvement on this imaging method was given by Chang and Edwards [240,241]; see Fig. 23. By delaying the read gradient echo a time tϕ relative to the spin echo, the chemical phase shift is encoded without the spatially varying background phase error. Accordingly, only Nϕ = 4 images are required to obtain accurate oil and water maps; specifically, ∆ωtϕ = −π /2, 0, π /2, and π , where the difference in chemical shift ∆ω = 2π (ν1 − ν2 ). This phase encoded chemical sensitive imaging technique is well suited to core analysis as good fluid-phase discrimination is achieved even for broad spectral lines that overlap. When spectroscopic resolution is not available, the easiest method of fluid-phase discrimination is to ensure only one fluid contains resonant spins. An obvious example is the displacement of brine by air, as imaged in capillary pressure measurements; see Section 4.4. If the rock contains two spin-bearing liquids, such as brine and oil, then the aqueous-phase may be replaced with D2 O as deuterium has a much lower gyromagnetic ratio than hydrogen and is therefore not observed in a typical 1 H measurement [242]. D2 O has also been used for miscible fluid displacement studies (i.e., D2 O displacing H2 O) [243], although care must be taken in such experiments to distinguish displacement from 1 H/2 H chemical exchange. A mixture of oil and heavy water does allow for the direct detection of 1 H (oil) and 2 H (brine) using a dual resonance rf probe [244]. Alternatively, the oil-phase could be fluorinated allowing signal to be detected from either 1 H or 19 F [173,216]. Oil displacement can also be observed through the 13 C nucleus [63], although the low natural abundance of the NMR sensitive carbon isotope makes detection difficult without isotopic enrichment or signal enhancement by polarization transfer [64]. Otherwise, the brine can be imaged independently by detection of the salt content through 23 Na [65,66] or 7 Li imaging [245]. Relaxation time is sensitive to chemistry and so T1 or T2 mapping can be used, under favorable circumstances, to identify the spatial distribution of oil and water. If one of the fluid-phases has a single, well-defined T1 relaxation time, this
J. Mitchell et al. / Physics Reports (
)
–
27
Fig. 23. Pulse sequence schematic for chemical selective imaging via phase encoding. An additional delay tϕ is introduced into a slice-selective spin-warp imaging sequence, ensuring the gradient echo (centered at time te + tϕ ) is delayed relative to the spin echo (centered at time te ). By varying tϕ in successive experiments, images of different chemical components are obtained.
component can be nulled so only the other fluid-phase is detected, as demonstrated in [246,247]; the T1 null is based on fat suppression in medical imaging [156]. Steady state free precession (SSFP) provides a similar suppression of a T2 component [248]. However, in rocks where T1 and T2 relaxation times are typically described by multi-exponential distributions, nulling of a single relaxation time component has limited applicability. More appropriately, spatially resolved T1 recovery or T2 decay curves can be fitted to obtain the intensity of the different relaxation time components. By attributing particular ranges of relaxation time to each fluid-phase, chemically specific maps are inferred [28,249]. The current approach is the generation of a relaxation time distribution for each pixel/voxel. The ranges of T1 or T2 are then integrated to provide improved intensities as described in [250]; examples of this method are presented in Section 5. When discrimination of the fluid-phases is not possible using the inherent relaxation time contrast, doping agents may be added to either the aqueous or organic fluidphase [251]. The use of doping agents also provides a route to fast imaging by reducing the T1 recovery time [252]. However, doping agents must be used with care in core analysis as the paramagnetic ions (typically, Mn2+ or Gd3+ ) will adsorb on a rock surface, and notably clays in sandstones [253], resulting in a time-varying concentration of doping agent and a loss of fluid-phase discrimination [254]. Chelating agents can be added to prevent adsorption, although these can chemically erode carbonate formations even under alkaline conditions [255]. When relaxation time contrast is not available, diffusion coefficients provide an alternative method of fluid-phase discrimination [256]. 4.3.1. Monitoring imbibition and recovery processes Spontaneous imbibition of liquids into rock occurs due to capillary action. Typically this process is slow (on the order of hours to days) and is therefore suitable for monitoring with MRI. The earliest 1D MRI of rocks addressed this question of capillary imbibition of water against air in a variety of materials [15]. Later, 2D imaging was used to monitor in situ uptake of water in limestone against gravity [257]. More recently, 1D MRI has been used to validate models of water imbibition into homogeneous rock subjected to pore-surface treatments [258]. Of more interest to the petroleum industry, the forced displacement of oil by water has been studied in 1D using various chemically selective techniques by Fordham et al. [239,259]. Similar 1D displacement experiments, imaged in 2D, are shown in [260]. Baldwin et al. [261] have demonstrated a 2D MRI equivalent of the standard Amott measurement of wettability index. The volume of oil displaced from a core-plug immersed in water is determined as a function of time; a limiting factor in this technique is oil viscosity, as high viscosity oils remain attached to the rock surface, even though they have been displaced from the pores. The MR images reveal the presence of surface oil left behind by displaced oil, providing an accurate measure of recovered oil volume. The reverse process of rock drying has also been studied [262]. Another imbibition process that has been studied by MRI is drilling mud invasion. Water or oil based muds are used down hole to facilitate cutting of the rock face and act as a lubricant and coolant for the drill bits. Drilling may be performed over-balanced with sufficient mud pressure to prevent kicks, although the drilling fluid will be forced into the surrounding formation under such conditions. Significant mud invasion prevents accurate assessment of oil saturation by logging tools. Similarly, laboratory measurements on extracted rock cores will be influenced by mud invasion. Limiting mud invasion whilst maintaining smooth cutting is an ongoing concern for well-log analysts. Fordham et al. used T1 weighting and profiling to monitor the ingress of clays into rock cores [16,17], whereas Straley et al. performed similar experiments with chemically sensitive imaging [263–265]. More recently, studies have been performed to assess the significance of mud invasion on MRI core analysis [266–268]. Forced displacement of oil by brine is a topic of considerable interest in core analysis as this mechanism underpins the secondary recovery process (water flood) from many reservoirs world-wide. Understanding of oil recovery processes at the laboratory-scale is critical for generating accurate models to predict field-scale performance in the reservoir. Some of the earliest MRI studies of forced oil displacement were conducted on a whole-body medical system operating at only ν0 = 10.6 MHz; this system was later modified to operate at ν0 = 21.3 MHz, thereby providing improved image
28
J. Mitchell et al. / Physics Reports (
)
–
Fig. 24. Maps of D2 O brine concentration calculated from 2D MRI images of remaining oil or H2 O brine under steady (Exp. 1) and unsteady (Exp. 4, 7) flow conditions. The formation of viscous instabilities is clear in the latter examples. Flow of solvent occurred from left (I, inlet) to right in all cases. In Exp. 1, the solvent injection front is distorted by gravity. Source: Images reproduced with permission from the Society of Petroleum Engineers [269].
quality of oil recovery from Berea sandstone [251,271,272], where paramagnetic doping of the aqueous fluid-phase was employed to provide chemical selectivity. Numerous displacement studies have been conducted using relaxation time weighting [18], CSI [273–276], or D2 O substitution [277–279] to distinguish the oil and aqueous fluid-phases. Determination of accurate flow parameters, notably relative permeability, is a topic of importance in core analysis as such data are used to define and validate numerical simulations of oil recovery. MRI has been employed in various ways to improve permeability estimates. Kulkarni et al. combined 1D profiles of remaining oil with pressure drop measurements to infer relative permeabilities [280,281]. Other authors have attempted similar calculations by imaging non-uniform saturation states [282–284]. Most forced displacement measurements are performed in 1D, although 2 and 3D images are useful for observing non-uniform displacement, as illustrated in Fig. 24 [269]. The possible existence of viscous instabilities, as observed in these MRI studies, has important implications for the injection of polymers [134,285,286] or gels [217] to improve oil recovery in EOR strategies. MRI studies of oil recovery from core-plugs are also used to validate numerical simulations [287–289]. Forced displacement of oil through fractured core-plugs has been a topic of some interest for MRI experimentalists, possibly because the fractures are readily visible in images or profiles [290–293]. Miscible liquid displacement in fractures has also been studied [294], as well as gas transport in chalk with styolites and a gouge-filled fracture [295]. Examples of oil recovery imaged in a fractured (composite) plug are shown in Fig. 25. Fluid in the fracture is identifiable by relaxation contrast, as the free-diffusing liquid will have considerably longer T1 and T2 relaxation times than liquid confined in the surrounding pore matrix [19,296]. MRI has been combined with nuclear tracer imaging (NTI) to monitor transport across a fracture [297,298]. Different transport properties are observed depending on whether the fracture is open or closed [299], and on wettability conditions [300,301]. The latter work was continued in [270,302] with improved image quality. Examples of direct gas imaging in rocks are scarce, although displacement of liquid by gas can be imaged readily, assuming the gas provides no significant NMR signal. Enhanced recovery of oil by supercritical CO2 has been studied [304,305], plus mixtures of supercritical CO2 and brine [306]. Salt (halite) precipitation after supercritical CO2 injection has been inferred from dry zones in a Berea sandstone [307]. Oil recovery through a fractured chalk was compared during gas (nitrogen) and water floods [308]; the oil recovery from the water flood was better than from the gas flood as the gas more easily by-passed the remaining oil. Gas hydrate formation and recovery is another area of interest. Gas hydrates (clathrates) consist of gas molecules (including light hydrocarbons, H2 S, and CO2 ) trapped in a frozen structure (cage) of water molecules. Hydrates form typically in deep marine and permafrost environments. Gas hydrates have a sufficiently short T2∗ not to be observed in most imaging techniques [309], whereas mobile water and methane are detectable. Consequently, gas hydrate formation is inferred from a loss of signal in MRI [303,310–312]. Example images following the formation of a gas hydrate in a Bentheim sandstone plug are shown in Fig. 26. 4.3.2. Flow mapping The concepts of phase-sensitive encoding for molecular diffusion and displacement discussed in Section 2.2 are also applied to measurements of velocity. When imaging flow through a porous medium such as rock, the displacement observed over the observation time is not scalable to velocity (as it is in pipe flow) due to the tortuosity of the pore network.
J. Mitchell et al. / Physics Reports (
)
–
29
Fig. 25. MRI of progressive oil recovery (A–F) through artificial rock fractures (composite plug). The gray scale ranges from black (no signal) to white (maximum signal). The main images show the transverse fracture and the inserts show the sagittal fracture. The aqueous fluid-phase (D2 O brine, no signal) was injected under constant pressure. The oil (decane) is clearly visible in the fracture (high signal intensity). Air bubbles become trapped during the course of the flood (example bubble circled in images B–F). Source: Images reproduced with permission from Springer [270].
Fig. 26. 3D MRI images of remaining free water and methane gas during hydrate formation in a Bentheim sandstone plug. Images acquired at experimental times of (a) 0 h, (b) 20 h, (c) 32 h and (d) 90 h. Color scale ranges from transparent (no signal) to dark blue (low signal intensity) to light blue/green (high signal intensity). No signal is obtained from the hydrate. MRI data were analyzed to estimate rate of hydrate formation. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Source: Images reproduced with permission from the Society of Petrophysicists and Well Log Analysts [303].
30
J. Mitchell et al. / Physics Reports (
)
–
Fig. 27. Pulse sequence schematic for displacement (velocity) encoding. The thin and thick lines indicate 90° and 180° rf pulses, respectively; the gray rectangle is a homospoil gradient. A series of 2Ng gradient pulses are applied — each of duration tδ and magnitude gv — interspersed between 180° rf refocusing pulses. The spins accumulate a phase shift proportional to their displacement. At the end of the sequence, the spins are stored on the z-axis by a 90° rf pulse. Any phase sensitive imaging sequence can be appended after the flow encoding; here, PEPI from Fig. 14 is suggested after [189].
Fig. 28. Displacement maps of water flowing in a Bentheimer core-plug obtained using flow encoded PEPI. Mean volumetric flow rates are indicated in each figure, with the flow rate increasing from (A–E). (F) and (G) are repeat measurements at constant flow rate, with the difference map given in (H). Variations in the difference map indicate that the flow behavior is not at steady-state. Source: Figure reproduced with permission from Elsevier [313].
Flow (displacement) encoding is achieved typically prior to the application of a standard imaging sequence. A displacement map can be obtained directly by virtue of a suitable phase pre-conditioning sequence designed to minimize signal loss due to diffusion in internal gradients; one version of this pre-conditioning sequence is shown in Fig. 27 [189,313]. A series of 2Ng gradient pulses, each of duration tδ and amplitude gv are applied between 180° rf refocusing pulses to encode for spin displacement. The total cycle (time between 90° excitation and z-storage rf pulses) generates a displacement sensitive phase accumulation of
ϕ=γ
t
gv t ′ v t ′ dt ′ ,
(35)
0
where t ′ = 2Ng tδ is the phase encoding time and v is the mean displacement over unit time (equivalent to velocity as t ′ → 0). The total phase shift due to flow is thus
ϕ = 4γ gv v Ng2 tδ (tδ + 2tδ1 ) ,
(36)
with tδ 1 being the time interval between a 180° rf refocusing pulse and the adjacent leading gradient edge (assuming the rf pulses are symmetric about the gradient lobes). The image generated following the velocity encoding will contain phase information that allows a mean velocity to be calculated for each voxel. Images obtained with flow encoded PEPI are shown in Fig. 28 with further examples in [313–317]. The flow encoding in Fig. 27 is equivalent to using the PGSE sequence in Fig. 6(a) with a long echo time. The reduced dephasing time achieved with multiple rf pulses in Fig. 27 provides some resilience to signal loss from displacement through internal gradients. The bipolar APGSTE sequence in Fig. 6(c) is more appropriate for flow encoding in porous media for the reasons discussed in Section 2.2. The net phase shift (proportional to displacement) is used to determine an apparent velocity based on the observation time t∆ . Accordingly, this technique has been used to provide spatially resolved displacement measurements (i.e., flow propagators) [318] and velocity maps [43] of liquid flowing in core-plugs. The gradient amplitude (max) is ramped between ±gv in sequential experiments, where the phase accumulated is
ϕ = γ gv t∆ tδ v,
(37)
and tδ is the gradient pulse duration as defined in Fig. 6. The differential of Eq. (37) provides dϕ dgv
= γ t∆ tδ v,
(38)
J. Mitchell et al. / Physics Reports (
)
–
31
and so the displacement per unit time (velocity) is obtained from the phase shift as a function of gradient strength. To obtain apparent velocity vectors, it is necessary to increment gv on all three axes (x, y, z) in sequential experiments. Edwards et al. [319] and Merrill [184,320] used this method to acquire displacement maps with spin-warp imaging. In the latter work, significant errors were found in the largest apparent mean velocities calculated from MRI, possibly due to spin relaxation or an inappropriate range of gradient amplitude. The flow displacements accessible by MRI are typically much larger than those encountered in a reservoir. Notwithstanding, flow imaging is useful in core analysis for determining the dominant transport pathways and connectivity in the rock formation [8,182], and for estimating local permeability [43]. 4.4. Capillary pressure measurements Determination of capillary pressure is important in core analysis as fluid transport and saturation states are affected by capillarity-driven processes [1]. Capillary pressure is generated by the contact of wetting and non-wetting fluids at a solid interface. ‘‘Capillary pressure curves’’ describe the variation in capillary pressure as a function of fluid saturation state and are used in core analysis to calculate petrophysical quantities such as irreducible water saturation, free water level, and residual oil saturation in the reservoir. Oil recovery predictions and simulations rely on capillary pressure measurements and therefore accurate determination of capillary pressure curves has significant industrial implications. A conventional capillary pressure curve is obtained using a centrifuge, where a core-plug is spun at a range of rotation speeds until the liquid distribution in the rock reaches equilibrium. The volume of liquid ejected (e.g., either oil or brine) is recorded as a function of centrifugal force [321,322]. To avoid demounting the plug to determine the volume of liquid displaced, a petrophysical centrifuge is modified with the addition of a camera and a stroboscopic light to photograph the free liquid captured in the sample holder (also modified with a viewing window calibrated volumetrically). A good quality capillary pressure curve will consist of data acquired over at least fifteen centrifuge speeds and is therefore time consuming (on the order of weeks) to complete. In addition to the complexity of modifying a centrifuge, the determination of capillary pressure curves still contains unresolved issues concerning the interpretation and assumptions in the data analysis [323]. In the original interpretation, capillary pressure Pc is defined as the difference between the pressures of the wetting Pwetting and non-wetting Pnon−wetting fluids in contact with a curved surface, such as the wall of a pore, as [321]: Pc = Pnon−wetting − Pwetting ,
(39)
and is related to the pore throat diameter dt of a cylindrical pore by the Washburn equation [324]: Pc = 4
σ cos (θ ) dt
,
(40)
where σ is the interfacial tension between the two fluids and θ is the contact angle. For a rock partially saturated with water, the rate of change of capillary pressure with height hc above the free water level is dPc dhc
= ∆ρ ag ,
(41)
where ∆ρ is the density difference between the wetting and non-wetting fluids and ag is the acceleration due to gravity. In the case of forced displacement in a centrifuge this equation becomes dPc dhc
= ∆ρ ac ,
(42)
where the centrifugal acceleration ac associated with angular velocity ϖc at a radius rc is given by ac = −ϖc2 rc . Therefore, assuming the capillary pressure at the outlet face of the core is zero, the capillary pressure as a function of radius is given by Pc (rc ) =
1 2
2 2 ∆ρϖc2 rc2 − rc1 .
(43)
However, in the traditional centrifuge method, a range of capillary pressures will exist within the length of a core, yet only the average volume of liquid displaced is observed. Hence the average saturation S¯ is measured, and is assumed to be equal to an integral across all capillary pressures, thus S¯ =
1
rc2 − rc1
0 Pc (rc1 )
S {Pc (rc )} dPc (rc ). −∆ρϖc2 rc
(44)
Eq. (44) is rearranged into the Hassler–Brunner integral for average saturation:
¯ c (rc1 ) = cos2 SP
Pc (rc1 )
C
2
0
S {Pc (rc )}
1 − {Pc (rc )/Pc (rc1 )} sin2 (C )
dPc (rc ),
(45)
32
J. Mitchell et al. / Physics Reports (
)
–
in which C = cos−1 (rc1 /rc2 ). Due to the use of average saturation, Eq. (45) is useful only when rc2 ≈ rc1 , i.e., in a very short core-plug. As Eqs. (41)–(45) are used to calculate Pc at each centrifuge speed, there is a recognized flaw in the method. Other approximate solutions to Eq. (45) have been suggested and are reviewed in [322,323], although these still contain assumptions and limitations. MRI profiling, combined with centrifugation, offers an improvement on the conventional centrifuge method by providing the spatial distribution of liquid saturation S (rc ) in a core-plug [38,176,275,325]. The profiling removes the necessity of attempting to deduce the capillary pressure from the average saturation via Eq. (45) and overcomes the requirement for the false assumption that rc1 ≈ rc2 . Each pixel in the saturation profile is small, thereby ensuring that the condition rc1 ≈ rc2 is met. As the MRI profile is acquired in the direction of the centrifugal acceleration, a considerable time-saving can be achieved: each pixel is associated with a different centrifugal acceleration according to Eq. (42) even though the centrifuge is operating only at a single speed. Therefore, it is possible to obtain the capillary pressure Pc (rc ) at all saturations S = 0 → 1 from a single profile. It is necessary to select the duration of centrifugation and centrifuge speed precisely to obtain such a profile that spans all saturation states. The experimental parameters will vary with geophysical rock properties (i.e., porosity and permeability) and fluid type. Calibrating the centrifuge speed/duration therefore requires a trial and error approach. To simplify the estimation of appropriate centrifugation parameters for a given sample, it is useful to convert the capillary pressure into the dimensionless Leverett J-function [326]: J (S ) =
√
κ/φ . σ cos (θ )
Pc (S )
(46)
The scaling law in Eq. (46) reflects the natural scaling of capillary pressure according to
σ cos (θ ) Pc (S ) = J (S ) √ ,
(47)
κ/φ
√
where (κ/φ) is a simple proxy for pore throat size based upon macroscopic, Darcy-scale parameters that are measured in routine core analysis. The J-function derives from the capillary bundle concept of a porous medium coupled with the Washburn equation (40). The value of J provides a convenient means of scaling the gross effects of interfacial tension, contact angle, and pore size, as a dimensionless capillary pressure characteristic. Hence porous media with geometrically similar pore matrices are expected to have identical J-functions. (max) reached at the In the context of a centrifuge capillary pressure determination, the maximum capillary pressure pc inner face of the core-plug is given by Eq. (43) with rc = rc1 . The scaled maximum is given by J
(max)
2 2 ϖc2 ∆ρ rc2 − rc1 κ = . 2σ cos (θ) φ
(48)
This provides a means of estimating the required centrifugation speed ϖc for a particular rock of given φ and κ , in a centrifuge of given dimensions, and with a liquid of known surface tension and contact angle against the rock minerals. A small set of J values will usually span the entire capillary pressure range of interest. The centrifuge speeds are determined by rearranging Eq. (48) to
ϖc2 = 2
J (max) σ cos (θ ) 2 ∆ρ rc2 − rc1
2
φ k
.
(49)
Therefore, in practice, saturation profiles are acquired at three or four different spinning speeds (Leverett J-functions) to determine a complete capillary pressure curve. Despite having to spin the core-plug at several J values, the combined MRI-centrifuge protocol offers a considerable time-saving over conventional capillary pressure curve measurements (experimental duration of days rather than weeks). To image the liquid distribution in the plug, it is necessary to remove the sample from the centrifuge and place it in the MRI magnet. During this transition, and for the duration of the imaging experiment, fluid redistribution can occur in the rock; this is a particular concern in high permeability samples where gravity overcomes capillarity [327]. Spinler, Baldwin, et al. prevented fluid redistribution when measuring oil (octadecane) and brine capillary pressure curves by freezing the oil fraction prior to removing the sample from the centrifuge [41,328,329]. The sample was therefore stable for the duration of the 2D spin-warp imaging experiments. Only the liquid water was observed in their images. In related work, capillary pressure curves have been determined at different wettability conditions [330] and from macroscopic regions of differing porosity [329]. The details of the imaging protocol can influence the quality of the measured capillary pressure curves as discussed in [331]. MRI has also been used to confirm sample orientation in X-ray measurements of capillary pressure [332]. An important application of the oil/brine capillary pressure measurement is the determination of wettability. Wettability refers to the tendency of a pore surface to adsorb either aqueous or organic species. In terms of reservoir engineering, rock cores can be water wet (imbibes water preferentially), oil wet (imbibes oil preferentially), or have mixed wettability whereby different surfaces within the same sample will be water or oil wet. Various indexes of wettability are defined in the literature, all based upon oil/brine capillary pressures. One such index is the industry-standard US Bureau of Mines (USBM) wettability index, Ux [333]. The USBM wettability index is the difference in the natural logarithms of the areas under
J. Mitchell et al. / Physics Reports (
)
–
33
Fig. 29. Schematic of ideal oil/brine capillary pressure curves, plotted against water saturation Sw . The areas under the positive and negative curves are given by A1 and A2 , respectively, and are used to determine the USBM wettability index according to Eq. (50).
Fig. 30. Capillary pressure curve measurement of brine draining from a Berea sandstone plug, obtained by SPRITE profiling. The raw profiles are shown in (a): prior to centrifugation with water saturation Sw = 100 s.u. everywhere (◦), and after centrifugation with an average water saturation S¯w = 46.3 s.u. (△). By taking the ratio of the partially saturated (△) to fully saturated (◦) profiles on a pixel-by-pixel basis, the water saturation Sw (z ) is obtained as a function of distance z from the plug inlet face (z = 1.8 cm). Here, z is aligned along the plug as defined in Fig. 11(a). A saturation state of Sw = 100 s.u. (free water level) is maintained at the outlet by a foot-bath. The distance z is converted to capillary pressure Pc by Eq. (43), thus providing the capillary pressure curve in (b). Typically, several centrifuge spinning speeds are used to ensure a uniform distribution of data over the full (accessible) range of Sw . Source: Images reproduced with permission from the American Institute of Physics [33].
a capillary pressure curve for a complete oil/water displacement/imbibition cycle, also equal to the work required for the non-wetting fluid-phase to displace the wetting fluid-phase. The USBM index is therefore given by
Ux = ln
A1 A2
,
(50)
where A1 and A2 are the areas under the positive and negative capillary pressure curves as defined in Fig. 29. Nørgaard and other authors demonstrated an improved method for determining oil/brine capillary pressure curves using 1D CSI [39,40,276,334]. As fluid transport occurs primarily in one direction along the plug, there is no requirement for multidimensional image acquisitions, and the chemical shift dimension provided simultaneous determination of the oil and brine saturations. More recently, SPRITE has been used as the profiling sequence of choice in combined MRI-centrifuge studies to provide quantitative capillary pressure curves [33,34]. The method has been demonstrated on a low-field magnet, B0 = 0.18 T (6.2 MHz), for improved quantification in core analysis [35,36,214,335]. The main limitation of the SPRITE protocol for capillary pressure measurements is the insensitivity to different chemical species. The lack of simultaneous oil and brine detection is not an issue for air/oil or air/brine capillary pressure measurements, but curtails application to oil/brine studies as H2 O must be replaced by D2 O. The use of heavy water modifies the result as the centrifugal force is dependent on the difference in fluid densities. Nevertheless, the rapid determination of accurate single fluid-phase capillary pressure curves is important in core analysis. An example of a 1D air/brine capillary pressure measurement using SPRITE is given in Fig. 30. Further examples of capillary pressure curves, determined using a low-field spectrometer, are given in Section 5.4. When a single liquid, such as brine, is drained from a core-plug, the capillary pressure is determined by the force required to drive the liquid through pore throats. By rearranging the Washburn equation (40) it is possible to convert Pc to a pore throat diameter dt by dt = 4
σ cos (θ ) Pc
.
(51)
Pore throat size distributions are obtained in conventional core analysis using mercury intrusion porosimetry (MIP) [336]. Green demonstrated that the pore throat size distributions obtained from SPRITE MRI-centrifuge data are consist with those obtained by MIP [37].
34
J. Mitchell et al. / Physics Reports (
)
–
Fig. 31. Spatially resolved T2 distributions obtained in the laboratory and displayed in the style of a well-log, after [338]. The T2 distributions were obtained from a limestone plug saturated with a light Middle East crude oil (left) before and (middle) after injection of brine. A T2 cut-off (vertical red line) is positioned to discriminate between the two fluid-phases. Here, the brine has been doped with manganese to reduce the relaxation time; hence the brine signal appears to the left of the oil signal in the T2 distributions (brine: blue, oil:black). The central portion of the plug (thicker lines on the T2 distributions where 1.4 cm < y < 2.5 cm) is considered to be free from capillary end effects. The spatial variation of porosity is shown on the right (solid line), along with the oil saturation before and after brine flood (dashed lines). The blue rectangle represents the volume of oil recovered from the central portion of the plug. The arrow indicates the direction of brine flow. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
5. Extension to very low field strengths In core analysis, very low magnetic fields are preferred (B0 ≈ 50 mT) to provide comparability to well-logging tools. Additionally, the use of low-field magnets limits the magnitude of internal gradients [9]; hence more robust T2 relaxation time measurements are achieved [10]. Based on the literature survey in Fig. 3 (Section 1), we see that MRI experiments conducted at such low magnetic field strengths have recently become popular. Imaging at B0 ∼ 50 mT has been considered impractical due to the poor SNR inherent from such systems, even though historically, early images were acquired at such magnetic field strengths [51]. Notwithstanding, very low-field magnets, designed predominantly for use in core analysis, have until recently been restricted to bulk measurements of core-plugs. Now, dramatic developments in NMR hardware have provided a much-needed boost to the capabilities of these very-low field systems; SNR has been improved by new rf amplifier technology and the progression to fully digital spectrometers, removing the need for an analogue multiplier. The SNR obtained from a standard core-plug of volume V ∼ 100 cm3 is sufficient to allow spatial resolution on the millimeter scale (at least in 1D) with a sensitivity equivalent to better than 1 s.u. In this section, we review and demonstrate the latest developments in very low-field MRI as applied to core analysis. We now use the notation that the plug is aligned along the y-axis, according to Fig. 11(b). Unless otherwise stated, all data in this section were acquired using an Oxford Instruments’ ‘‘Big-2’’ magnet (B0 = 50 mT) controlled by a Maran DRX HF spectrometer, equivalent to the Geospec2 system shown in Fig. 2 (Section 1). 5.1. Relaxation time mapping One of the most useful profiling techniques applied at very low-field is a T2 map spatially resolved on the y-axis. The y–T2 map is a powerful measurement in core analysis as it provides data akin to those obtained from well-logs; example plots are given in Section 5.2. A NMR well-log displays a series of T2 distributions plotted as a function of reservoir height h, i.e., a h–T2 map. A core-plug profile obtained in the (y, log10 T2 ) plane under laboratory conditions is similarly a series of T2 distributions plotted as a function of plug height (here, y is the vertical dimension). An example of such a series of spatially resolved T2 distributions, obtained in the laboratory but displayed in a manner consistent with a well-log, is shown in Fig. 31. In the field, well-logs may be obtained before and after a flood to determine the change in fluid saturation in the vicinity of the well bore. Consistent data are obtained from a core-plug by forced displacement of oil by brine; comparing
J. Mitchell et al. / Physics Reports (
)
–
35
Table 1 Radio frequency phase cycle for the basic CPMG experiment, also suitable for the T2 -mapping pulse sequences illustrated in Fig. 32. 90° rf pulse
180° rf pulse
Acquisition
x x x¯ x¯ y y y¯ y¯
y y¯ y y¯ x x¯ x x¯
x x x¯ x¯ y y y¯ y¯
the results obtained at the two length-scales (well bore and core-plug) provides further information on the performance of the oil recovery process [337]. Hence, these data emphasize the potential synergy of laboratory-scale core analysis and field-scale well-logging. The data present in a y–T2 map can be simplified for further analysis or comparison by taking the projection (sum) of the data along each axis. A projection on the y axis (summed along the T2 dimension) will be a profile where the intensity may be calibrated in porosity units. A projection on the T2 axis (summed along the y direction) will be a so-called ‘‘stacked’’ T2 plot. This data format is immediately interpretable by log analysts because of the similarity to the conventional well-log display. The only significant difference between the y–T2 map and a well-log is the spatial resolution: in the laboratory, ∆y is on the order of a millimeter, whereas in the reservoir, ∆h is on the order of tens of centimeters. Current research is focused on rapid methods of obtaining y–T2 maps that are suitable for time-critical studies of liquid transport in rocks, such as in situ studies of oil recovery. Here, we present a comparison of typical results obtained using the frequency encoded multi-echo sequence, the original SESPI sequence, and the modified SESPI sequence at very lowfield to complement recent publications. The sequence of rf pulses in these y–T2 acquisitions is identical to that in a CPMG measurement; appropriate phase cycles for the rf pulses and acquisition are given in Table 1. Many T2 weighted imaging sequences have been presented in the medical literature, often incorporating complicated phase cycles, composite rf pulses, or homospoil gradients. Here we demonstrate only the basic multi-echo sequence, as presented in [28], and SEPSI sequences, as presented in [228], to show that acceptable results are achieved at low-field without such modifications. Therefore, y–T2 pulse sequences can be incorporated readily into multi-dimensional acquisitions where the last dimension is almost always acquired using a CPMG echo train, such as T1 –T2 or D –T2 correlations [59]. When utilizing the frequency encoded multi-echo sequence, the spatial dimension is incorporated into these multi-dimensional acquisitions without a significant increase in experimental time compared to conventional bulk measurements. As we have seen, the use of frequency encoding in petrophysical studies has been considered inappropriate due to local gradient distortions (including internal gradients) or the requirement to measure very short T2 components [228]. When using very low field magnets, we argue that frequency encoded profiles are suitable for many core-plug measurements because:
• the low spatial resolution means the acquisition window is narrow, ensuring T2∗ contrast is minimized; • the minimum echo time is equivalent to that readily achievable in the original SESPI experiment; • the T2∗ decay of a liquid saturated rock tends to be mono-exponential [213] allowing straight-forward correction of the echo shape in k-space;
• the use of low field magnets ensures that T2 is maximized by minimizing internal gradients [10]; • the high temporal resolution makes the measurement ideal for dynamic systems such as monitoring forced displacement of oil. The multi-echo and SESPI pulse sequences are shown in Fig. 32, the schematics having been adjusted to reflect the implementation on very low-field magnets. Trapezoidal gradient pulses are shown explicitly (it is good practice in all MRI experiments to use ramped gradients to minimize eddy current stabilization delays) to emphasize the proximity to the rf pulses. Careful positioning of the rf and gradient events allows the echo time to be minimized in each experiment. Processing of the y–T2 data is divided into the following steps. First, the k-space data are Fourier transformed. The FT is performed individually on the k-space data acquired at each of the n echo times. As such, the resultant data contain n profiles, the amplitude of which decreases exponentially with T2 . In the case of the frequency encoded k-space data, it is desirable to include a first order phase shift. This is achieved by barrel shifting the k-space echo such that the maximum signal intensity corresponds to k = 0 cm−1 . The SESPI experiment includes the acquisition of a datum at k = 0 cm−1 , preventing the need for this pre-processing step. Once the FT is complete, a zeroth order phase shift is applied to rotate the majority of the complex signal into the real channel. The imaginary channel can be used to estimate the SNR but is not required otherwise. An appropriate numerical inversion method is then applied to invert the entire y–T2 map in a single operation using the method for solving 2D Fredholm integral equations given in Section 2.4.1. The relevant 2D Fredholm
36
J. Mitchell et al. / Physics Reports (
)
–
Fig. 32. T2 mapping pulse sequences for very low-field core analysis applications. Three alternative pulse sequences are shown: (a) the frequency encoded multi-echo method (after Fig. 18), (b) the original SESPI method, and (c) the modified SESPI method (both after Fig. 21). The thin and thick vertical bars represent 90° and 180o rf pulses, respectively. In (a), the ‘‘read’’ gradient has an area equivalent to gtr . In (b) and (c), the ‘‘phase’’ gradient lobes have areas that range over ±g (max) tp . In all cases, n spin echoes are observed; in (a) the full echo shape is recorded, whereas in (b, c) only the echo amplitude is measured. These sequences are modified relative to those shown in Figs. 18 and 21 to minimize the echo time and other delays.
integral equation is b (nte , k) b(0, 0)
∞
∞
K0 (T2 , nte , y, k) f (T2 , y) d (log10 T2 ) d (y) + ε,
= −∞
(52)
−∞
where the kernel matrix K0 is separable under the Krönecker product K0 = K2 ⊗ K1 such that B = K1 FKT2 + E.
(53)
Here, K1 describes the T2 decay – being identical in functional form to Eq. (4) – and K2 describes the FT. Since the FT inverse problem is well-posed, it is more efficient to invert the spatial dimension separately using a fast Fourier transform (FFT) algorithm. When implementing the 2D inversion, K2 is then replaced by the identity matrix of equal size. This 2D inversion is preferable to inverting each 1D T2 decay separately as it provides consistent smoothing over the entire distribution and reduces noise-induced discontinuities between adjacent T2 distributions. 5.2. Comparison of data To compare the results obtained from the three imaging sequences in Fig. 32, an imaging phantom was constructed using three tubes of deionized water doped to varying degrees with a paramagnetic species to give a range of T2 contrast. The phantom components were deionized water (T2 = 2 s), NiCl2 solution (T2 = 0.16 s), and CuSO4 solution (T2 = 0.04 s). A schematic of the phantom is shown in Fig. 33. A small quantity of surfactant solution (soap) was added to the deionized water to reduce the curvature of the meniscus. The y–T2 maps obtained from the imaging phantom using the three alternative pulse sequences in Fig. 32 are shown in Fig. 34. In all cases ν0 = 2 MHz (1 H) and the recycle delay was tRD = 15 s. The frequency encoded y–T2 data were acquired with te = 2 ms, n = 4000 echoes, g = 2.74 G cm−1 , m = 64 points per echo with a dwell time of ∆t = 10 µs, a digital filter with bandwidth fBW ∼ 100 kHz, and 32 repeat scans. The total experiment time was 12 min. The SESPI y–T2 data (original and modified sequences) were acquired with tp = 400 µs, m = 64 phase encode increments spaced linearly between ±g (max) = ±2.22 G cm−1 , 11 points per echo with a dwell time of ∆t = 1 µs, a digital filter with bandwidth
J. Mitchell et al. / Physics Reports (
)
–
37
Fig. 33. Imaging phantom constructed from cylindrical tubes of (a) deionized water, (b) NiCl2 solution, and (c) CuSO4 solution. The y-axis is vertical. The total phantom height was y ≈ 4 cm.
fBW ∼ 2 kHz, and 4 repeat scans. In the original SESPI pulse sequence te = 2 ms and n = 4000 echoes; in the modified SESPI pulse sequence the initial echo time was te = 2 ms, the reduced echo time was te′ = 1 ms, and n = 8000 echoes. For both phase encoded sequences the total experiment time was 98 min. Despite having a reduced number of scans, the SESPI experiments have a significantly longer acquisition time than the frequency encoded equivalent, cf. 98 min versus 12 min. The frequency encoding y–T2 map, Fig. 34(a), features a T2 distribution with narrow peaks. The quality of the T2 decay is improved by the imaging gradients that act as homospoils to eliminate non-ideal spin behavior. The T2 distributions obtained using the original SESPI method, Fig. 34(b), and modified SESPI method, Fig. 34(c), are less well defined. This loss in T2 quality is particularly noticeable in the modified SESPI sequence due to the complicated spin dynamics. Returning to the frequency encoded y–T2 map, Fig. 34(a), we note that the T2 of the deionized water is reduced from T2 = 2 s to 1.8 s due to diffusion in the presence of the imaging gradients. Notwithstanding, this pulse sequence is still applicable to petrophysical samples as the tortuous pore structure of rock limits the free diffusion path length and hence the magnitude of this effect. It is expected that surface relaxation will be a more significant contributor to T2 than diffusion in the direction of the read gradient in such samples. If the diffusive contribution to T2 is considered important in a particular experiment, then a diffusive decay term can be included in the inversion process to obtain a corrected T2 distribution [125]. The frequency encoded y profile in Fig. 34(a) is distorted slightly due to the convolution with the spectral linewidth (T2∗ ). In core analysis applications, any such distortion would be deconvolved. For the purposes of comparing the data obtained from the three pulse sequences considered, we applied the same processing steps in the generation of each y–T2 map and therefore retain the pulse sequence-dependent artifacts in the data shown. Additionally, the T2∗ decay of this doped water phantom was not mono-exponential; consequently, the complicated deconvolution process was not attempted here. Nevertheless, the integral area under the profile is still quantitative even without this correction, as demonstrated in Section 5.3 and [134]. The image profile shape is less well defined in the original SESPI y–T2 map, Fig. 34(b), with distortions arising from imperfect matching of the positive and negative gradient lobes. The modified SESPI method provides the best profile shape, Fig. 34(c), as it relies on the application of only a single phase gradient pulse. However, the T2 dimension is over-smoothed by the inversion algorithm to accommodate the complicated spin dynamics (resulting in non-exponential decays) inherent in the pulse sequence design. An estimate of the SNR was calculated based on the T2 decay at, or close to, the origin of k-space. This can be achieved readily in the phase encoded SESPI sequences as a CPMG decay is acquired with g = 0 G cm−1 . However, in the frequency encoded data, the maximum amplitude signal acquired at the center of the echo has to be taken as an approximation for k = 0 cm−1 . On a per scan basis, the frequency encoded data were obtained with SNR = 30 and the phase encoded data were obtained with SNR = 180. This difference is due to the variation in digital filter bandwidth. Therefore, the frequency encoding method requires 36× more repeat scans than the phase encode methods to attain a given SNR. As this number of repetitions is less than the number of phase encode increments (m = 64), frequency encoding is faster. It is worth noting at this point that the absolute SNR becomes less significant in these measurements due to the data processing. The spin echoes may be window averaged in each CPMG-like decay to improve the effective SNR [92] and the numerical inversion used to
38
J. Mitchell et al. / Physics Reports (
)
–
Fig. 34. T2 maps acquired from the imaging phantom in Fig. 33 using various methods: (a) frequency encoding, (b) original SESPI, and (c) modified SESPI.
obtain the T2 dimension performs well in the presence of noise, as long as the noise is white and Gaussian with a mean of zero. Good quality y–T2 maps are still obtained even when the SNR is reduced, as demonstrated in Fig. 34. Additionally, we draw attention to the use of a custom narrowband digital filter in the SESPI measurement to maximize the SNR. The difference in SNR per scan between the frequency and phase encoded profiles would be much less if standard MRI filters were used for both acquisitions. When selecting the optimum pulse sequence for a specific application, a trade-off must be made between the quality of the T2 distribution and the profile shape. In core analysis studies, fluid-phase discrimination is based on a T2 cut-off and it is therefore imperative to obtain an accurate relaxation time measurement. As such, it is preferable to use frequency encoding when quantitative estimates of oil and brine volumes are required. In other applications, such as capillary pressure
J. Mitchell et al. / Physics Reports (
)
–
39
measurements, an undistorted profile shape is necessary, in which case the modified SESPI sequence would become the y–T2 encoding method of choice. As discussed earlier, if SNR is unimportant then the frequency encoded y–T2 data are acquired more rapidly than the phase encoded data. However, SESPI based y–T2 maps, under-sampled in k-space and reconstructed using CS (see Section 6.1), may provide an advantage in experimental time when a particular SNR is required, especially in low porosity rocks. 5.3. Monitoring forced displacement of oil As mentioned earlier in Section 4.3.1, T2 maps allow simultaneous quantification of the volume of oil and aqueous fluidphases in a rock, and are therefore valuable when monitoring forced displacement of oil [134,339]. The spatial dimension provides a visualization of the displacement mechanism that can be used to provide input parameters for reservoir-scale simulations of oil recovery. We demonstrate the use of y–T2 mapping to monitor oil displacement by brine injection in a model system of dodecane in sandstone. The experiment mimics the log-inject-log protocol of a piloting operation used to determine the efficacy of EOR agents in the reservoir [151]. Example y–T2 maps obtained using the frequency encoded multi-echo sequence in Fig. 32(a) during an oil recovery experiment are shown in Fig. 35. A water-wet core-plug of Bentheimer sandstone with dimensions 60 mm × 38 mm (length × diameter) was saturated with dodecane. The Bentheimer sandstone, favored previously in NMR studies due to its low paramagnetic ion content, has a porosity φ ≈ 0.24 and a gas permeability κ ≈ 3000 mD [10]. The total pore volume of the core-plug was approximately φ V = 15 cm3 . Known volumes of 2 wt% KCl brine were injected into the core-plug and y–T2 maps were acquired between these brine flood stages. The brine was doped with MnEDTA to provide T2 contrast relative to the dodecane. The y–T2 map of the dodecane saturated plug (initial state) is shown in Fig. 35(a) and is seen to be saturated uniformly with a monomodal T2 distribution. Flow occurred vertically from the inlet at y = −3 cm to the outlet at y = 3 cm. After 3 cm3 of brine had been injected, Fig. 35(b), the brine (short T2 component) is seen predominantly at the base of the core-plug (inlet) and the brine saturation decreases towards the top of the core-plug (outlet) as expected. The dodecane saturation shows the opposite trend and hence, as seen from the profile projection, the core-plug remains at a uniform saturation state overall. A T2 cut-off was located at T2 ≈ 350 ms to separate signal from the doped brine and dodecane. Once 15 cm3 of brine had been injected, Fig. 35(c), both the brine and dodecane are seen to be distributed uniformly throughout the core-plug. This suggests the brine has fingered through the dodecance but failed to displace all of the oil. The SNR decreased as the flood progressed because the electrically conductive brine provided a route for interference to reach the rf probe. This is an inherent problem in these type of flow-through experiments and can be addressed only by careful earthing of the flow lines [115]. The optimum T2 cut-off shifted slightly throughout the experiment; in Fig. 35(c) it is located at T2 ≈ 0.43 s due to a narrowing in the T2 component attributed to the oil. The y–T2 maps obtained during the core flood were analyzed to provide quantitative volumes of oil and brine present in the sandstone. A calibration signal was determined for known volumes of brine and oil, i.e., magnetization (arbitrary units) per unit volume of liquid [134]. Once this calibration was obtained, both the oil and brine volumes present in the core-plug were calculated for each stage of the flood. Using conventional gravimetric or volumetric techniques, this simultaneous determination of both fluid-phase saturation states can be achieved only by measuring the volume of oil and brine in the effluent and accounting accurately for all dead volumes in the flow lines. Effluent collection is complicated when studying heavy oils where precipitates may become lodged in the flow line, or live oils with a high volatile content. An additional disadvantage of a volumetric/gravimetric measurement to determine the remaining oil is the accumulation of errors in the estimation of the oil saturation. The y–T2 map provides an independent determination of volume for each fluid-phase at each stage of the flood. As the NMR measurement probes the saturating liquid directly, this measurement overcomes any difficulties in capturing and assaying the effluent in a quantitative manner. Overall, NMR estimates of remaining oil can be more accurate than determinations based on effluent analysis. A comparison of NMR and volumetric oil and brine saturations are shown in Fig. 36 where excellent agreement is seen between the two sets of data. Additionally, the NMR data confirmed (within experimental error) that the total saturation of the core-plug did not vary during the experiment. The NMR data provide a direct, simple, and robust measure of the oil and brine saturations compared to effluent analysis, and the NMR analysis does not require the assumption that the total saturation remains constant. Although the oil recovery studies demonstrated here and in [134,339] were conducted in a stop–start manner, the temporal resolution of the y–T2 mapping experiments is sufficient to allow continuous monitoring of saturation state during the displacement process (i.e., while fluids are flowing), particularly when using flow rates consistent with those encountered in oil fields some distance from the well bore. Typical interstitial flow velocities relevant to oil recovery are often on the order of v = 10−4 cm s−1 , allowing time for the acquisition of y–T2 maps without deleterious temporal blurring. SESPI has been used to monitor the recovery of a light Middle East crude oil from a microporous limestone core-plug during a continuous brine flood [338]. The cylindrical limestone plug had a porosity φ = 0.33 and a Klinkenberg gas permeability κ = 10−2 µm2 , with dimensions 37.5 mm × 50.9 mm (diameter × length); 1 PV ≡ 18.3 cm3 . The initial (initial) oil saturation (determined gravimetrically) was So = 76 s.u. The core-plug was held in a NMR-compatible core-holder at a temperature of 58 °C and an isostatic confining pressure of 2 MPa for the duration of the flood. Brine was injected at a volumetric flow rate of V˙ = 1.4 × 10−3 cm3 s−1 , equivalent to an interstitial velocity of v = 1 ft day−1 in the formation. SESPI y–T2 maps were acquired continually throughout the flood with m = 32 phase gradient increments, a maximum
40
J. Mitchell et al. / Physics Reports (
)
–
Fig. 35. Frequency encoded y–T2 maps of oil recovery by brine flood from a water-wet Bentheimer sandstone core-plug saturated initially with dodecane. Brine (2 wt% KCl solution) doped with MnEDTA was injected into the core-plug. The y–T2 maps were acquired following injection of (a) 0 cm3 , (b) 3 cm3 , and (c) 15 cm3 of doped brine. Flow occurred vertically from the core-plug inlet face at y = −3 cm to the outlet face at y = 3 cm. The vertical line in (b, c) indicates the optimum T2 cut-off for discrimination between oil (long T2 ) and brine (short T2 ) fluid-phases in the analysis.
gradient strength g (max) = 3.6 G cm−1 , and gradient pulses of duration tp = 300 µs. The profiles had FOV = 6.5 cm and resolution ∆y ≈ 2 mm. The echo time was te = 2 ms, n = 128 echoes were observed, and 16 repeat scans were co-added to improve the SNR. Each y–T2 map was acquired in 45 min. The brine was estimated to move approximately 7 mm through the plug during one SESPI acquisition, introducing some blurring into the profiles. A T2 cut-off was defined at T2 = 0.03 s between the brine doped with Mn2aq+ (short T2 ) and the oil (long T2 ). The time evolution of the spatially resolved oil and brine saturations are shown in Fig. 37(a) and (b), respectively. The total liquid saturation is shown in Fig. 37(c), where it is seen to remain at S = 100 s.u. throughout the flood. The remaining oil at the end of the flood was distributed non-uniformly through the plug, with So increasing towards the outlet face. The high oil saturation observed at the outlet was attributed
J. Mitchell et al. / Physics Reports (
)
–
41
Fig. 36. Quantitative oil recovery monitoring by NMR. Oil and brine saturations determined from the y–T2 maps in Fig. 35 are compared against a volumetric determination of oil recovered from a dodecane saturated Bentheimer sandstone core-plug. The symbols correspond to NMR data (squares) and volumetric data (triangles) for the oil (filled symbols) and brine (open symbols). The total liquid volume (gray squares) was obtained from the NMR data. The volumetric determination of brine saturation was calculated from the recovered oil fraction, assuming the core-plug was saturated fully at all times. Error bars (equivalent to ±0.5 cm3 per measurement) are shown for the volumetric data only; errors on the NMR data are comparable to the marker size.
to a capillary end effect. In conventional bulk studies of saturation (determined either by bulk NMR analysis or gravimetric assays of recovered oil) such variation in oil saturation cannot be assessed. Continuous crude oil recovery has also been monitored in combination with chemical EOR agents [337]. An alkali surfactant (AS) preparation was chosen to improve the oil recovery. In these experiments, microporous limestone plugs were again used, here with typical porosity φ = 0.29, Klinkenberg gas permeability κ = 6 × 10−3 µm2 , and 1 PV ≡ 15.2 cm3 . The plugs had average dimensions 38 mm × 50 mm (diameter × length). Heterogeneity maps of these plugs are shown in Fig. 17, obtained at intermediate-field. Each plug was mounted in the NMR compatible core-holder at elevated temperature and pressure for the duration of the flood, which consisted of injecting 15 PV brine, followed by a variable volume of AS solution, then additional brine to a total flood volume of 25 PV. All fluids were injected at a flow rate equivalent to an interstitial velocity of v = 1 ft day−1 . Due to the presence of the alkali, the brine could not be doped with Mn2aq+ (precipitates at high pH) and a chelating agent such as EDTA would chemically erode the limestone even under alkali conditions [255]. Therefore, D2 O was used in all aqueous fluid-phases to suppress the 1 H signal. T2 maps of the oil were obtained using the frequency encoded multi-echo method with te = 1000 µs, n = 1024 echoes, g = 2.7 G cm−1 , m = 64 data points, and a dwell time ∆t = 10 µs. The profiles had FOV = 8.6 cm and resolution ∆y ≈ 1.3 mm. Each y–T2 map was acquired in 15 min so the injected fluids were estimated to move approximately 3 mm through the plug during one multi-echo acquisition. No significant changes in the T2 of the oil were noticed during the flood; measurement of the T2 dimension was still important as it allowed a quantitative determination of the oil saturation without T2 weighting. The time evolution of the spatially resolved oil saturations are shown in Fig. 38 for injected AS volumes of (a) 8 PV and (b) 2.24 PV. In both cases, a high oil saturation remains at the outlet during the initial brine flood stage (volume pumped <15.5 PV). On introducing the AS solution, additional oil is liberated from the pores and forced through the rock; the process leading to this localized oil resaturation and then desaturation of oil is referred to as the transport of an oil bank. As the AS solution reaches the end of the plug, it sweeps out the oil trapped at the outlet by capillarity. Further injection of AS solution, Fig. 38(a), results in a continued recovery of a small volume of oil, notably at the inlet face. When studying short core-plugs, end effects due to geometry, viscous instabilities, or capillarity can accommodate a significant portion of the sample volume. The spatial resolution achievable with MRI is therefore essential to elucidate the behavior of the fluids in a section of the plug unperturbed by such end effects. As demonstrated in [134,337], the spatial resolution allows the portion of the plug unaffected by end effects to be studied independently, thereby providing an indication of fluid behavior in a much larger (infinite) volume; ideally, this volume will be representative of the reservoirscale. 5.4. Capillary pressure measurements The significance of combining MRI with centrifuge capillary pressure measurements was discussed in Section 4.4. Currently, quantitative SPRITE profiling at low to intermediate field strengths is the preference for this type of profiling. However, the inclusion of a T2 dimension allows the centrifugal drainage of fluid from pores of different size to be observed. Implementation on a very low-field magnet ensures the technique can be ported readily into existing NMR core analysis work flows [340]. An example of a single fluid-phase (brine) forced drainage experiment monitored by y–T2 mapping is shown in Fig. 39 for Leverett-J values of (a) J = 0, (b) 2, and (c) 8. The rock sample was an Estaillades limestone core-plug of dimensions 50 mm × 38 mm (length × diameter). Estaillades limestone has a porosity of φ ≈ 0.3 and a permeability κ ≈ 100 µm2 .
42
J. Mitchell et al. / Physics Reports (
)
–
Fig. 37. Time evolution of spatially resolved (a) crude oil, (b) brine, and (c) total fluid-phase saturations in a microporous limestone plug, after [338]. T2 maps were acquired continuously using SESPI during the oil recovery process and a T2 cut-off was used to provide oil/brine fluid-phase discrimination. The inlet face of the plug was located at y = 0 cm; the brine flow was from left to right. Bulk fluid saturations are projected on the right for clarity. A capillary end effect (high oil saturation) is observed near the plug outlet. The typical error in the saturations (color scale) is less than ±5 s.u..
This rock has a bimodal pore size distribution, evident in the T2 distribution shown in Fig. 39(a) and confirmed by MIP. The T2 dimension in Fig. 39 has been rescaled to pore body diameter db ≡ 6ρ2 T2 assuming the industry standard limestone surface relaxivity of ρ2 = 1 µm s−1 . The y–T2 maps reveal the different spatial distributions of brine obtained in the micropores (db < 1 µm) and macropores (db > 1 µm) after centrifugation at three different spinning speeds. As the core-plug drains, it is clear from the T2 distributions that the large pores empty first, leaving only the microporosity saturated at J = 2. The y–T2 map in Fig. 39(c) shows that the loss of the long T2 component is consistent along the length of the core-plug. This behavior would not be evident from the profiles alone whereas the y–T2 maps provide a straight-forward interpretation. The analysis of the y–T2 maps suggests that liquid transport occurs through the well-connected macroporosity independent of the microporosity. This is in agreement with previous studies of Estaillades limestone where no diffusive exchange of water was observed between the micropores and macropores on the time-scale of an NMR experiment [149]. When the profile projections in Fig. 39 are analyzed using the method described in [33] (see Fig. 30 for a visual demonstration), then a capillary pressure curve for the entire rock core-plug is obtained, see Fig. 40(a). The T2 dimension in Fig. 39 allows the drainage of the micropores and macropores to be considered independently. It is seen from Fig. 40(a) that the maximum spinning speed was insufficient to drain the micropores. Therefore we determine a capillary pressure curve for the macroporosity only, see Fig. 40(b). Note that this macroporosity capillary pressure curve is not just a simple re-scaling of the bulk measurement in Fig. 40(a). The ability to consider the capillary pressure response of the macropores independently is important as these large pores dominate liquid transport, and hence oil recovery, in these bimodal pore
J. Mitchell et al. / Physics Reports (
)
–
43
Fig. 38. Time evolution of spatially resolved crude oil saturations in microporous limestone plugs during an EOR process, after [337]. T2 maps were acquired continuously using frequency encoded multi-echo profiling during the flood. The inlet face of the plugs was located at y = 0 cm in each case; fluid flow was from left to right. AS injection occurred between the horizontal (white) dashed lines. The variation in bulk oil saturation is shown projected on the right of each plot. For both plugs, an initial piston-like displacement of oil is observed followed by a gradual recovery during the brine flood; elevated oil saturation is observed at the outlet: this is a capillary end effect. On AS injection, increased mobility of the oil results in the transport of an oil bank through the plug; the AS also removes the capillary end effect. Note that the temporal resolution of these experiments is a factor 3 better than in the equivalent SESPI acquisitions (brine only flood) shown in Fig. 37. The typical error in the saturations (color scale) is less than ±5 s.u.
networks. If a larger range of spinning speeds were available, separate capillary pressure curves would be obtained for both the macropores and micropores. It is known from Eq. (40) that capillary pressure Pc is scalable to pore throat diameter, dt . Likewise, T2 is scalable to pore body diameter db from Eq. (65). Assuming σ = 72 mN m−1 is the air/brine interfacial tension, θ = 60° is the air/brine contact angle on the limestone, and ρ2 = 1 µm s−1 , then the average pore throat diameter ⟨dt ⟩ is correlated with the average pore body size ⟨db ⟩ to obtain a body-to-throat ratio (BTR) of the macropores, see Fig. 40(c). A positive correlation is observed, suggesting that the largest throats are associated with the largest pores and vice versa, as expected [82]. In the smallest macropores, BTR > 10, whereas in the largest macropores, BTR → 1. Determining the BTR is important for understanding well-logs in microporous limestone reservoirs [341], and the T2 mapping MRI-centrifuge protocol allows a direct comparison of pore body and pore throat size from a single measurement. We suggest that continued investigation into y–T2 maps of forced imbibition or drainage of liquid will provide new insights into transport mechanisms through rock. In the case of the MRI-centrifuge experiment, the T2 dimension will enable simultaneous determination of oil/brine capillary pressure curves by virtue of a T2 cut-off, as demonstrated in Section 5.3. 5.5. Multi-dimensional relaxation and diffusion mapping Correlations of NMR parameters in 2D have become an important resource in well-logging and core analysis. For example, diffusion-relaxation D –T2 correlations provide improved oil/brine discrimination where 1D T2 distributions alone are insufficient [134]. Rauschhuber and Hirasaki demonstrated a spatially resolved D –T2 measurement of core-plugs at low-field by replacing the CPMG portion of the pulse sequence with the frequency encoded multi-echo sequence [59]. The multi-echo profiling sequence can be substituted for the CPMG portion in any of the D –T2 pulse sequences in Fig. 9 or 10. Relaxation correlations are also useful, particularly in the laboratory where T1 –T2 plots are used to determine the ratio T1 /T2 , known to be related to the strength of interaction between the imbibed liquid and the rock mineral surface [142–144]. Previously, spatially resolved T2 measurements were used as a local probe of wettability index at very low-field [342]. However, T2 is inherently sensitive to pore geometry and this must be taken into account when scaling to some (arbitrary) wettability index. T1 /T2 has the advantage of removing (to a good approximation) the geometry term [142], thus providing a more direct link to wettability. Here we demonstrate the spatially resolved T1 –T2 measurement. The T1 dimension is encoded using the standard saturation recovery method [67], see Fig. 4(b), with τ1 varied logarithmically in 16 steps from τ1 = 0.1 to 10 s. A recycle delay of tRD = 1 s was used. The multi-echo profiling portion of the sequence was acquired with te = 1000 µs, n = 1024 echoes, g = 2.7 G cm−1 , m = 64 data points, and a dwell time ∆t = 10 µs. The total acquisition time of the y–T1 –T2 map was 1 h. The y–T1 –T2 experiment is demonstrated here using a composite plug of Estaillades limestone formed from two short plugs (length 2.5 cm) stacked end-to-end. One half of the composite plug was saturated with 2 wt% KCl brine and the other half was saturated with hexadecane (oil). The clean limestone is water-wet and so a stronger surface interaction is expected for the brine than the oil. Example T1 –T2 correlations are shown in Fig. 41 typical of (a) the brine saturated half of the composite plug and (b) the oil saturated half. In the water saturated portion, Fig. 41(a), an average (brine) T1 /T2 = 2.0 is observed, indicating a strong interaction between the brine and the carbonate surface. The individual T1 and T2 distributions are bimodal and reflect the pore size distribution of the rock. In the oil saturated portion, Fig. 41(b), an (oil) average T1 /T2 ≃ 1.0 due to the weak interaction with the surface. The individual T1 and T2 distributions are monomodal and do not contain information on the pore structure, perhaps because the oil does not penetrate the micropores. It is clear from Fig. 41(c) that a bulk T1 –T2 correlation of the entire composite plug does not provide sufficient oil and water fluid-phase discrimination.
44
J. Mitchell et al. / Physics Reports (
)
–
Fig. 39. T2 maps of forced drainage of a brine saturated Estaillades limestone core-plug. The core-plug was spun in a standard centrifuge at speeds corresponding to Leverett-J numbers of (a) J = 0, (b) 2, and (c) 8. The T2 dimension is ‘‘unit-normalized’’ assuming ρ = 1 µm s−1 .
A spatially resolved indicator of wettability will be of great benefit in core analysis. Sandstones typically exhibit mixed wetting properties and the ability to assess the heterogeneity associated with the liquid/solid interactions will provide additional insight into oil recovery processes. It is common practice in core analysis to attempt to restore the reservoir wettability conditions after Soxhlet cleaning of the plugs. The reliability of wettability restoration methods is an ongoing topic of debate within the industry so experiments that can determine accurately the degree of wettability alteration are of interest to petrophysicists. Furthermore, some chemical EOR agents (e.g., surfactants) function by switching the wettability of the mineral surface from oil-wet to water-wet. A spatially resolved T1 /T2 profile would allow the efficacy of such agents to (brine) be examined directly. An example T1 /T2 profile is shown in Fig. 41(d) for the composite limestone plug with T1 /T2 ∼2 (oil)
(where −2.5 cm < y < 0 cm) and T1 /T2 ∼ 1 (where 0 cm < y < 2.5 cm). A clear distinction is visible between the brine saturated and oil saturated portions and shows overall that the rock surface is strongly water-wet. For these data, acquired at ν0 = 2 MHz, the full range of sensitivity available in T1 /T2 between a ‘‘wetting’’ and ‘‘non-wetting’’ liquid is a factor 2. Improved sensitivity can be achieved by acquiring the T1 /T2 profile at higher field strengths, e.g. ν0 = 10 MHz, as T1 /T2 increases with Larmor frequency [142]. Small errors in the T2 dimension due to diffusion in the (weak) internal gradients generated at such modest field strengths are removed by determining the diffusive contribution to the CPMG decay [125].
J. Mitchell et al. / Physics Reports (
)
–
45
Fig. 40. Capillary pressure curves for brine drainage from an Estaillades limestone core-plug. Utilizing the y–T2 maps shown in Fig. 39, we obtain (a) a capillary pressure curve for the entire plug, and (b) a capillary pressure curve for the macropores only. Smoothed capillary pressure curves (solid lines) are shown as a guide to the eye. A BTR plot is shown in (c) correlating dt (from Pc ) against db (from T2 ) for the macropores.
6. Advanced imaging techniques Depending on the precise application of MRI in core analysis, obtaining a high resolution 3D image of a core-plug can be both time consuming and unnecessary. In this section we consider recent advances in MRI technique that offer alternatives to regular k-space sampling; these advanced techniques allow us to make measurements that were hitherto considered impractical using MR. We will consider two examples: compressed sensing MRI and Bayesian MR. In the former, image acquisition is accelerated by under-sampling in k-space so images are acquired with a higher temporal resolution or with a higher spatial resolution in a given acquisition time. In Bayesian MR, a different approach entirely is used. Bayesian MR is useful when, for example, details of a particular structural length-scale within a system are required. Bayesian MR identifies the relevant information directly in the k-space data without need for acquiring a complete image. We illustrate this method with an application to grain sizing in core-plugs. 6.1. Compressed sensing Compressed sensing (CS) is a technique for obtaining an image from an incomplete set of k-space data, acquired using a specific (pseudo-random) sampling scheme. If we have prior knowledge that the image can be represented sparsely in some transform domain, e.g., following the discrete cosine, wavelet, or curvelet transform, then we can use the theory of CS to reconstruct the image from a subsample of the full data. If the data (in this case k-space points) are under-sampled incoherently (meaning the sampling scheme is not sparse in the transform domain [343]) then the image is obtained using an appropriate non-linear reconstruction method. The CS reconstruction replaces the FT stage in image processing.
46
J. Mitchell et al. / Physics Reports (
)
–
Fig. 41. Spatially resolved T1 –T2 correlations of a composite Estaillades limestone plug containing 2 wt% KCl brine (−2.5 cm < y < 0 cm) and hexadecane (0 cm < y < 2.5 cm). Typical T1 –T2 correlations are shown for (a) the brine saturated portion and (b) the oil saturated portion. A standard bulk T1 –T2 plot is shown in (c) for comparison. In (a–c) the dashed diagonal line indicates T1 /T2 = 1 and the solid diagonal line indicates T1 /T2 = 2. The contour intervals are the same in both plots. A T1 /T2 profile extracted from the y–T1 –T2 map is shown in (d).
CS reconstruction provides a significant improvement to phase encoded MRI acquisitions in particular, due to the nature of the sampling scheme (one wavenumber per measurement) where any k-space position can be readily selected. SPIbased acquisitions provide an obvious example of an imaging method wherein k-space wavenumbers may be sampled incoherently [222,344]; at the end of this section we demonstrate the combination of CS reconstruction and T2 -mapping with SESPI to improve the temporal resolution of the measurement at very low-field. Elsewhere, CS has been applied to other MR measurements, such as multi-dimensional relaxation time [345] and flow propagator acquisitions [346]. We now describe a CS reconstruction method after [60,61], derived from basis pursuit [347,348]. This reconstruction uses the ℓ1 -norm to characterize the sparsity of the reconstructed image. First the nD k-space data are stacked into a vector b. The image to be reconstructed is likewise stacked as a vector f and this is transformed into the sparse domain by operator Ψ . The k-space data vector b is obtained from the image by the under-sampled Fourier transform F . The reconstruction is achieved by solving the constrained optimization problem min ∥Ψ f∥1 , s.t. ∥F f − b∥2 ≤ ε,
(54)
where ε is a threshold that is set equal to the expected noise level and the ℓ1 -norm, ∥f∥1 = j fj , is a surrogate for sparsity. Therefore, minimizing the objective in Eq. (54) provides an image that is consistent with the acquired k-space data but has the most sparse representation achievable in the transform domain. Eq. (54) is appealing because it is a convex optimization problem and therefore can be readily solved using a number of algorithms. Eq. (54) can either be solved directly or in its unconstrained Lagrangian form:
arg min ∥F f − b∥2 + λ1 ∥Ψ f∥1 f
(55)
using methods including projected conjugate gradients algorithm [348], iterative thresholding [349], interior-point methods [350], and Bregman iterations [351]. The basic CS algorithm can be improved by including additional prior knowledge. One example often used is the combination of the spatial finite-differences (total variation) transform along with another sparse transform [352], such as a wavelet transform. The optimization problem then becomes min (∥Ψa f∥1 + λ2 ∥Ψb f∥1 ) , s.t. ∥F f − b∥2 ≤ ε,
(56)
J. Mitchell et al. / Physics Reports (
)
–
47
Fig. 42. Demonstration of 1D compressed sensing applied to the original SESPI method: (a) under-sampled k-space data (real component), equal to 50% of the full acquisition; (b) the CS reconstructed y–T2 map is a good approximation to the equivalent map, sampled fully, shown in Fig. 34(b).
where Ψa is the wavelet operator, Ψb is the spatial finite-differences operator, and the constant λ2 trades sparsity in each of the two domains. As with Eq. (54), Eq. (56) may be solved in an unconstrained Lagrangian form: arg min ∥F f − b∥2 + λa ∥Ψa f∥1 + λb ∥Ψb f∥1 . f
(57)
Alternatively, if the sample position within the image is known, this information can be incorporated into the reconstruction. The image is masked by a binary stacked vector M (with elements having a value 1 where sample is known to exist and zero otherwise), so Eq. (56) has an additional constraint s.t. ∥(1 − M) ◦ f∥2 ≤ ε,
(58)
where the ◦ operator indicates an element-by-element vector multiplication. The unconstrained Lagrangian form is: arg min ∥F f − b∥2 + λa ∥Ψa f∥1 + λb ∥Ψb f∥1 + λc ∥(1 − M) ◦ f∥2 , f
(59)
with the relative weighting in each term controlled by λa , λb , and λc . When the required result is known to be real, as with many MRI acquisitions, it is simpler to reconstruct only the real part of a phase corrected signal. However, if complex data are required, as in phase sensitive velocity imaging, the real and imaginary components are reconstructed separately and the phase is calculated afterwards. Methods for reconstructing phase sensitive k-space data are described in [61]. We now present a demonstration of CS reconstruction applied to low-field MRI acquisitions. Specifically, we consider the use of CS to improve the temporal resolution of the SESPI measurement. To reduce the acquisition time of the SESPI measurement, CS can be used to reconstruct y–T2 data under-sampled in the spatial dimension. To demonstrate this process, a repeat y–T2 map was acquired using the original SESPI sequence with 50% of k-space selected incoherently. The sampled points were chosen from a Gaussian probability distribution to ensure the center of k-space is sampled completely. The under-sampled 1D k-space data are shown in Fig. 42(a). A combination of total variation and Haar wavelet transforms were used according to Eq. (57). The total experiment time was reduced to 49 min (from 98 min) and as such the under-sampled phase encoded data would be faster to acquire than the frequency encoded data if a particular SNR was required. In 2D or 3D imaging, it is typical to sample less than 30% of k-space and still recover the image without noticeable distortion [61]. However, in the low resolution 1D data acquired here, a 30% sampling scheme resulted in significant artifacts in the reconstructed profile. A minimum of 50% sampling was found to be necessary to avoid these reconstruction artifacts. The CS reconstructed y–T2 map in Fig. 42(b) is a good representation of the fully sampled map shown in Fig. 34(b), see Section 5.2. The peaks in the T2 distribution are broadened slightly following the CS reconstruction and some minor oscillations are smoothed out of the profile as a consequence of the wavelet transform. Overall there is no significant difference between the fully sampled and under-sampled y–T2 maps. CS reconstruction is an important addition to the MRI toolbox whenever phase encoded acquisitions are necessary and can be used either to provide a significant time saving or improve SNR without increasing the total experiment duration. 6.2. Grain sizing: a Bayesian-MR approach Conventional grain sizing in core analysis is achieved using scanning electron microscopy (SEM) or X-ray microcomputed tomography (XMT). Both of these methods examine small sample volumes (V ∼ 10 mm3 ) that will not be representative of a rock formation with macroscopic (millimeter-scale) structural heterogeneity. Example images reconstructed from XMT data are shown in Fig. 43. Ideally, one would like to obtain 3D MR images on core-plugs with
48
J. Mitchell et al. / Physics Reports (
)
–
Fig. 43. Example images reconstructed from XMT data are shown for (a) Bentheimer sandstone and (b) Ketton limestone. The Bayesian MR measurement is sensitive to the large oolithic structures in the limestone. Images previously shown in [62].
sufficient spatial resolution to identify individual grains. However, obtaining such images from a liquid-saturated rock is prohibited by the very long acquisition times. For example, a spin-warp experiment required to obtain a 3D resolution ∆ {x, y, z } ∼ 10 µm with the required SNR would have a total duration in excess of 50 years. The fast imaging techniques described in Section 4 provide lower SNR and the effective resolution is reduced (blurred) by the relaxation weighting in these multi-echo sequences. Furthermore, the sample volume would be severely limited to satisfy the Fourier relationship between FOV and voxel size, and so an entire core-plug could not be imaged. Recently, a Bayesian inference approach was introduced to overcome these limitations by extracting a grain size distribution directly from 1D k-space data [62]. Mansfield and Grannell’s original work on MRI drew analogies between k-space and scattering patterns as observed in X-ray diffraction experiments [168]. Although the practicalities of imaging with magnetic resonance and X-rays are very different, k-space data do contain information on structural length-scales and diffraction patterns have been observed in MRI studies of regular, repeating granular materials. The Bayesian inference approach described in [62] demonstrates that structural information is also extracted from measurements of disordered granular materials. Hence, Bayesian inference can be used to predict a likely distribution of grain sizes that gives rise to the interference pattern observed in k-space. Bayesian inference was combined previously with MR to extract spectra from noisy data [353], flow velocities from under-sampled data [354], and bubble size distributions in a bubble column with a gas-voidage higher than can be explored by conventional optical techniques [355]. Bayes’ theorem [356] can be used to relate the probability of a hypothesis to a set of measured data. Formally, this is expressed by the posterior probability distribution P(ϑ|ι), which describes the probability of a hypothesis ϑ given the measurements ι. The posterior probability distribution is calculated from P(ϑ|ι) =
P(ι|ϑ)P(ϑ) P(ι)
,
(60)
where P(ι|ϑ) is the likelihood function describing the probability of observing the data given that the hypothesis ϑ is true and P(ϑ) is the prior, which describes our knowledge of the system before the measurements were made. The denominator in Eq. (60), P(ι), is called the evidence and describes the probability of measuring any set of data ι. In applications of parameter estimation, the evidence is a constant that acts only to normalize the probability distribution, and can be ignored. Eq. (60) is used to estimate the grain size distribution by parameterizing the grain size distribution in terms of, for example, a normal distribution. In this case, the data are the k-space acquisitions, so ι ≡ {|b(k)|}, and the hypotheses to consider are the unique parameters of the grain size distribution, so ϑ ≡ {¯r , ςr }, where r¯ is the mean grain size radius and ςr is the standard deviation of the distribution. The Bayesian-MR approach to grain sizing is limited only by magnet geometry and gradient strength, and can therefore be applied to standard sized plugs (V ∼ 106 mm3 ). As noted in Section 3.2.1, the range of wavenumbers spanned in k-space is given by k = γ gt /2π . The MR signal obtained from a liquid saturated core-plug is bs (k) = b0 (k) − b(k), where b0 (k) is the signal obtained from a bulk liquid sample of equal dimension to the core-plug, and b(k) corresponds here to the interference pattern due to the rock grains. The range of k to be analyzed is ideally chosen such that b0 (k) → 0 whilst |b(k)| > 0, so that |bs (k)| ≈ |b(k)| for all k of interest. As the Bayesian analysis is applied to the magnitude of the k-space data, it is sufficient to assume that |bs (k)| is equivalent to the interference pattern if an appropriate range of k can be selected; otherwise, it is necessary to account for b0 (k) in the data processing [355]. To determine the most likely grain size distribution that gives rise to the observed interference pattern, it is necessary to model the k-space data. Based on information extracted from XMT images, Fig. 43, the grain size distributions of the rocks were expected to be normal distributions and therefore modeled by
(r − r¯ )2 P (r ; r¯ , ςr ) = √ exp − . 2ςr2 2π ςr 1
(61)
For the Bayesian analysis of grain size, Eq. (60) becomes P (¯r , ςr ||b(k)| ) ∝ P ( |b(k)|| r¯ , ςr ) P (¯r ) P (ςr ) ,
(62)
J. Mitchell et al. / Physics Reports (
)
–
49
Fig. 44. (a) Typical k-space data (solid line) acquired for a frequency-encoded 1D profile of a water saturated Bentheimer sandstone. The noise floor estimate is shown by the dashed line. The fluctuations in the signal correspond to phase interference and are greater than the noise for k ≤ 104 m−1 . (b) Grain size distribution from Bayesian analysis of k-space data of Bentheimer sandstone (blue line) and Ketton limestone (red line), compared to XMT data (symbols) of the same. Data previously shown in [62]. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
where P (¯r ) and P (ςr ) are uniform, independent priors in the limit that ςr ≤ r¯ /2. The likelihood function relating the magnitude of the k-space data to the parameters of the grain size distribution is the Rayleigh distribution [357] P ( |b(k)|| r¯ , ςr ) ≡ P ( |b(k)|| R (k; r¯ , ςr )) =
|b(k)| − |b(k)|2 exp . R 2 (k; r¯ , ςr ) 2R 2 (k; r¯ , ςr )
(63)
The Rayleigh (width) parameter R is determined from the expected (average) form of the interference pattern such that R 2 (k; r¯ , ςr ) =
1 2
L/2
N
P (r ; r¯ , ςr ) |p (r ; k)|2 dr ,
(64)
0
where p (r ; k) is the simulated k-space response from a single grain and N is the number of grains in the sample of length L. A validation of this choice of likelihood function is given in [355]. Using Eq. (63), the posterior probability is calculated for all combinations of r¯ and ςr , and the optimum grain size distribution parameters are estimated from the mean of the posterior distribution. Holland et al. demonstrated the Bayesian grain sizing technique on water-saturated Bentheimer sandstone and Ketton limestone core-plugs at B0 = 1 T [62]. The method should be equally applicable on low-field and very low-field instrumentation, within the limits of achievable SNR. Appropriately high k-space wavenumbers are sampled to encompass the length-scales being measured. An example of typical k-space data and the resulting grain size distributions for the two rocks are given in Fig. 44(a) and (b), respectively. The grain size distribution can be used to predict the corresponding pore body size distribution following a MonteCarlo approach [358]. By simulating 2D close-packed clusters of three grains with random size (chosen based on the grain radius probability distribution), a characteristic length-scale (equivalent to pore body diameter) db is determined for the interstice between the circular projection of the spheres. The interstice of area A is approximated by an equilateral triangle, √ so db ≈ 2 A/2. Arrangements of more than three grains could be simulated when appropriate [359]. The process for extracting a pore size distribution is illustrated in Fig. 45. We know from Eq. (5) that the observed T2 relaxation time is considered proportional to pore body size as [4] 1 T2
=
1 T2,bulk
+ ρ2
6 db
.
(65)
Determining the surface relaxivity term ρ2 for sandstones is notoriously difficult as gas (krypton) adsorption measurements provide unreliable surface area results in the presence of clays. Accordingly, Holland et al. used the pore size distribution inferred from the Bayesian-MR grain size measurement to rescale the T2 distribution obtained for water saturated Bentheimer sandstone to pore size. They obtained a surface relaxivity of ρ2 = 8 µm s−1 [62]; not unreasonable for such a rock. Scaling of carbonate T2 distributions to pore size can be achieved reliably using surface area measurements from gas adsorption. Allen et al. explored the surface relaxivities of a selection of carbonate formations and discovered that although the value of ρ2 varied, it was always small and typically in the range ρ2 = 1–3 µm s−1 [82]. A value of ρ2 = 1 µm s−1 is considered the industry standard for carbonate formations when detailed analysis is unavailable. Holland et al. used this value to rescale the T2 distribution for water saturated Ketton limestone to pore size, in order to compare against the pore size distribution derived from the Bayesian-MR grain size distribution. The T2 distribution obtained on a B0 = 50 mT magnet was seen to be bimodal as expected; however, due to the weak scaling of ρ2 in Eq. (65), the long T2 component was dominated by the bulk 1/T2,bulk relaxation rate and was therefore not sensitive to the large pores between the ooids.
50
J. Mitchell et al. / Physics Reports (
)
–
Fig. 45. Generation of pore size distribution based on Bayesian-MR grain size distribution, demonstrated for Bentheimer sandstone: (a) the grain size distribution (repeated from Fig. 44) is used to generate (b) a close packed arrangement of three spherical grains projected in 2D. The interstice (blue) between the grains is used to estimate a characteristic pore body diameter db . This process is repeated 106 times to create the representative pore size distribution in (c, points). A T2 distribution (c, line), obtained using a very low-field magnet, is scaled to match the pore size distribution with ρ2 = 8 µm s−1 . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 46. Pore body size distribution obtained from a combination of T2 relaxation time and Bayesian-MR grain size measurements for Ketton limestone: (a) the T2 distribution is rescaled to pore diameter db using the industry standard ρ2 = 1 µm s−1 . The short T2 component reflects the microporosity (intragranular porosity) whereas the long T2 component is only sensitive to the bulk liquid relaxation time and not the pore size. The Bayesian-MR derived pore size distribution (b) is sensitive only to the macroporosity (intergranular porosity). Combining the microporosity distribution from T2 and the macroporosity distribution from Bayesian-MR, the entire pore body size distribution is obtained in (c).
The short T2 component was sensitive to the microporosity within individual ooids that could not be probed by the BayesianMR methodology. Therefore, by combining the pore size distributions obtained from T2 and Bayesian-MR, the entire pore size distribution of the rock was elucidated, as illustrated in Fig. 46. At present, Bayesian-MR grain sizing of rocks has been demonstrated only on an intermediate field magnet, although it is expected to perform well on low-field systems. The technique is robust in the sense that it is reasonably insensitive to the liquid in the pores; grain size can be determined even from a mix of oil and brine, as long as the proton density is spatially uniform. 7. Conclusions Throughout this review, the power of magnetic resonance imaging as a tool for core analysis has been irrefutably demonstrated, whether as a qualitative indicator of structural heterogeneity or a quantitative measure of fluid saturation and transport. The ability to spatially resolve the distribution of fluid-phases (oil, brine, gas) on timescales of seconds to minutes
J. Mitchell et al. / Physics Reports (
)
–
51
enables in situ studies of oil recovery processes. Accordingly, laboratory-scale MRI is used to determine the efficacy of EOR agents – e.g., chemical (surfactants, polymers), supercritical CO2 , or miscible gas – prior to piloting operations in a reservoir. The ability to probe nuclei other than 1 H allows studies of miscible liquid injection, such as monitoring 23 Na to study oil recovery with modified salinity brines. MRI has also been combined with conventional core analysis techniques, such as centrifugation, to rapidly determine capillary pressure curves. Ongoing research into spatially resolved relaxation times and diffusion coefficients is enabling the generation of wettability maps: another important tool for core analysis, leading to improved understanding of oil recovery processes. Recent advances in MRI techniques will enable the determination of rock structure (grain and pore size distributions) and fluid transport (flow propagators) on the same core-plug using a combination of MRI experiments. These laboratory-scale measurements are critical for improving models to predict field-scale performance of oil recovery strategies, and it is envisioned that the results will provide future insights into the structure-transport relationships that govern fluid flow in reservoirs. The early days of MRI in core analysis were dominated by intermediate-field measurements using commercial medical spectrometers. Unfortunately, the inherent magnetic susceptibility contrast in rocks (more significant in iron- or manganese-bearing sandstone formations) leads to enhanced relaxation time weighting in images as the magnetic field strength increases. Intermediate-field images, while informative in terms of structural heterogeneity on the millimeter scale, are not quantitative when fast (medical) MRI sequences are employed. Quantification of liquid saturations is achieved at intermediate to high field strengths using specialized MRI acquisitions such as SPRITE. These modern techniques are often limited, however, by hardware constraints, spatial resolution (image blurring), and a lack of inherent fluid-phase sensitivity. SPRITE is well suited to monitoring single fluid-phases, such as brine saturation when measuring an air/brine capillary pressure curve. It is less appropriate for simultaneous monitoring of oil and brine in a forced displacement (oil recovery) experiment. Recent developments in NMR hardware means that quantitative image profiles can be acquired on very low-field magnets, designed to provide consistent measurements between laboratory-scale core analysis and ‘‘reservoir-scale’’ well-logging. Relaxation time or diffusion contrast can be incorporated into profiling sequences to provide fluid-phase discrimination. A y–T2 map, for example, can be acquired on a B0 = 50 mT magnet (equivalent to ν0 = 2 MHz for 1 H) in a few minutes, allowing continuous monitoring of oil recovery from a core-plug at realistic reservoir flow rates. The inclusion of a T2 dimension is also significant in capillary pressure measurements as a body-to-throat ratio is determined from the simultaneous (correlated) measure of pore body size (from T2 ) and pore throat size (from Pc ). Low-field MRI has the additional advantage of minimizing magnetic susceptibility contrast, allowing quantitative fluid saturations to be determined using conventional (rapid) imaging techniques. As spectrometer hardware continues to evolve, it is envisioned that more sophisticated pulse sequences and techniques will be ported from the medical community to core analysis laboratories. High resolution MRI is already achievable at modest field strengths of B0 ∼ 200 mT, offering the potential for 3D monitoring of oil recovery processes, and hence providing further insight into structure-transport relationships in these complicated materials. The continued evolution of MRI technology, while of great significance to the core analysis community, will also impact other areas. For example, rapid high resolution imaging at low magnetic fields will expand the possibilities for medical diagnosis using safe, low cost permanent magnet technology. Acknowledgments The authors thank Andy Sederman and Alex Taylor for helpful discussions. For assistance with sample preparation and experiments, the authors thank John Staniland, Laurence Hawkes, and Alex Wilson. For financial support, J.M. thanks Schlumberger Gould Research; T.C.C. thanks the EPSRC (grant EP/F041772/1); D.J.H. thanks Microsoft Research Connections and the EPSRC (grant EP/F047991/1). Appendix A. Nomenclature
ag ; ac A
α
b; b; B B0 B1 ∆B0 ∇ B0
β c c
Acceleration due to gravity; centrifuge Area Small tip angle rf pulse Data function/scalar; vector; matrix Magnitude of static magnetic field Magnitude of applied (rf) magnetic field Variation in static magnetic field Magnetic field gradient Power-law scalar Scalar, defines inversion (c = 2) or saturation (c = 1) of spins Vector of fitting parameters (continued on next page)
52
C dt db D
δ ε; e f ; f; F
ˆf; Fˆ F
φ ϕ
g; g g0 geff gv
γ h hc
θ ϑ i
ι j J k; k ∆k K; K
κ ℓs ℓg ℓe
L L
λ m M n nc Nϕ Ng Nv
ν ∆ν p Pc P
Ψ
q; q r rc r R R
ρ ∆ρ
S Sp /Vp
σ
J. Mitchell et al. / Physics Reports (
)
–
Constant Pore throat diameter Pore body diameter Diffusion coefficient Chemical shift Experimental error (noise) function/scalar; vector Distribution function; vector; matrix Fitted distribution vector; matrix Fourier transform operator Porosity Phase angle Magnetic field gradient magnitude; vector Background magnetic field gradient magnitude Effective internal magnetic field gradient magnitude Magnetic field gradient magnitude for velocity encoding Gyromagnetic ratio Height in reservoir Height above free water Contact angle Hypothesis (Bayes’ √ theorem) Imaginary unit ( −1) Experimental data (Bayes’ theorem) Index counter Leverett-J function Wavenumber magnitude; vector Resolution in k-space Kernel function; matrix Permeability Characteristic pore size Dephasing length scale Diffusion path length over time te Sample length Smoothing operator Weighting parameter Number of k-space data Image mask vector Loop counter, specifically number of echoes in CPMG Number of saturation rf pulse loops Number of phase increments Number of gradient pulses for velocity encoding Number of pairs of data used to generate variogram Resonant frequency: Frequency bandwidth of slice Projection of grains in k-space Capillary pressure Probability (distribution) Sparse domain transform operator Magnetization wave number; vector Radius Radius of rotation in centrifuge Position vector Displacement vector Rayleigh parameter Surface relaxivity Fluid density difference Saturation state Surface-to-volume ratio Interfacial tension (continued on next page)
J. Mitchell et al. / Physics Reports (
ς t; t tc tδ t δ 1 ,2 t∆ te texp tϕ tp tr tRD ts tse tstore T1 T2 T2∗ T2,eff
τse τ1 τ2 ∆t
u; u Ux
Υ v V V˙
ω ∆ω W0
ϖc
x, y, z ∆x, ∆y, ∆z ∆ xs
ξ Ξ
)
–
53
Standard deviation Experimental time scalar; vector Constant time Gradient pulse duration Delay before/after gradient pulse Observation time Echo time in CPMG Total experiment (scan) time Phase evolution time Phase gradient duration Read gradient duration Recycle delay between scans Slice gradient duration Spin echo time, specifically in diffusion experiments Storage time of spins on z-axis Longitudinal relaxation time Transverse relaxation time Transverse relaxation time in an inhomogeneous magnetic field Effective transverse relaxation time Half spin echo time (diffusion and imaging) Longitudinal recovery time Half echo time (CPMG) Dwell time Fitting parameter scalar; vector USBM wettability index Variogram function Interstitial velocity Total core-plug volume Volumetric flow rate Larmor frequency Chemical shift difference Precision matrix Angular velocity of centrifuge Cartesian coordinates Spatial resolution Slice thickness in image space Pixel position Variogram parameter.
Appendix B. Abbreviations
CPMG CS CSI FOV FT ILT MRI MSME NMR PEPI p.u. PV RARE
Carr–Purcell–Meiboom–Gill (pulse sequence) Compressed sensing (reconstruction) Chemical shift imaging Field of view Fourier transform Inverse Laplace transform Magnetic resonance imaging Multi-slice–multi-echo (pulse sequence) Nuclear magnetic resonance π -echo planar imaging (pulse sequence) Porosity units Pore volume Rapid acquisition with relaxation enhancement (pulse sequence) (continued on next page)
54
J. Mitchell et al. / Physics Reports (
)
–
rf SESPI SPI SPRITE
Radio frequency Spin echo single point imaging (pulse sequence) Single point imaging (pulse sequence) Single point ramped imaging with T1 enhancement (pulse sequence) STRAFI Stray-field imaging s.u. Saturation units USBM United States Bureau of Mines.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40]
F.A.L. Dullien, Porous Media: Fluid Transport and Pore Structure, Academic, New York, 1991. R.L. Kleinberg, J.A. Jackson, An introduction to the history of NMR well logging, Concept. Magn. Reson. 13 (2001) 340–342. R.L. Kleinberg, Well logging overview, Concept. Magn. Reson. 13 (2001) 342–343. D.P. Gallegos, K. Munn, D.M. Smith, D.L. Stermer, A NMR technique for the analysis of pore structure — application to materials with well-defined pore structure, J. Colloid Interf. Sci. 119 (1987) 127–140. P. Callaghan, Principles of Nuclear Magnetic Resonance Microscopy, Clarendon, Oxford, 1991. F. Dalitz, M. Cudaj, M. Maiwald, G. Guthausen, Process and reaction monitoring by low-field NMR spectroscopy, Prog. Nucl. Magn. Reson. Spect. 60 (2012) 52–70. L.D. Brown, NMR imaging: Principles and recent progress, SCA Conference Paper Number 12, 1988. J.W. Gleeson, D.E. Woessner, C.F. Jordan Jr., NMR imaging of pore structures in limestones, SPE Form. Eval. (1993) 123–127. M.D. Hürlimann, Effective gradients in porous media due to susceptibility differences, J. Magn. Reson. 131 (1998) 232–240. J. Mitchell, T.C. Chandrasekera, M.L. Johns, L.F. Gladden, E.J. Fordham, Nuclear magnetic resonance relaxation and diffusion in the presence of internal gradients: the effect of magnetic field strength, Phys. Rev. E 81 (2010) 026101. I. Foley, S.A. Farooqui, R.L. Kleinberg, Effect of paramagnetic ions on NMR relaxation of fluids at solid surfaces, J. Magn. Reson. Ser. A 123 (1996) 95–104. C.E. Renshaw, J.S. Dadakis, S.R. Brown, Measuring fracture apertures: a comparison of methods, Geophys. Res. Lett. 27 (2000) 289–292. J. Arnold, C. Clauser, R. Pechnig, S. Anferova, V. Anferova, B. Blümich, Porosity and permeability from mobile NMR core-scanning, Petrophysics 47 (2006) 306–314. S. Anferova, V. Anferova, J. Arnold, E. Talnishnikh, M.A. Voda, K. Kupferschläger, B. Blümler, C. Clauser, B. Blümich, Improved Halbach sensor for NMR scanning of drill cores, Magn. Reson. Imaging 25 (2007) 474–480. R.J. Gummerson, C. Hall, W.D. Hoff, R. Hawkes, G.N. Holland, W.S. Moore, Unsaturated water flow in porous materials observed by NMR imaging, Nature 281 (1979) 56–57. E.J. Fordham, M.A. Horsfield, L.D. Hall, G.C. Maitland, Depth filtration of clay in rock cores observed by one-dimensional 1 H NMR imaging, J. Colloid Interf. Sci. 156 (1993) 253–255. E.J. Fordham, A. Sezginer, L.D. Hall, Imaging multiexponential relaxation in the (y, loge T1 ) plane, with application to clay filtration in rock cores, J. Magn. Reson. A 113 (1995) 139–150. G.C. Borgia, A. Brancolini, A. Camanzi, G. Maddinelli, Capillary water determination in core plugs: a combined study based on imaging techniques and relaxation analysis, Magn. Reson. Imaging 12 (1994) 221–224. G.C. Borgia, V. Bortolotti, P. Fantazzini, Magnetic resonance relaxation-tomography to assess fractures induced in vugular carbonate cores, SPE, 56787, 1999, ATC, Houston, Texas, 3–6 October, 1999. K.E. Washburn, C.D. Eccles, P.T. Callaghan, The dependence on magnetic field strength of correlated internal gradient relaxation time distributions in heterogeneous materials, J. Magn. Reson. 194 (2008) 33–40. A. Matteson, J.P. Tomanic, M.M. Herron, D.F. Allen, W.E. Kenyon, NMR relaxation of clay/brine mixtures, SPE Reservoir Eval. & Eng. 3 (2000) 408–413. S.V. Dvinskikh, K. Szutkowski, I. Furó, MRI profiles over very wide concentration ranges: application to swelling of a bentonite clay, J. Magn. Reson. 198 (2009) 146–150. A.J. Fagan, N. Nestle, D.J. Lurie, Continuous wave MRI diffusion study of water in bentonite clay, Magn. Reson. Imaging 23 (2005) 317–319. D.N. Guilfoyle, P. Mansfield, Rapid volumetric NMR imaging of fluids in porous solids using a 3D π -EPI (PEPI) hybrid, J. Magn. Reson. Ser. A 119 (1996) 151–156. A.M. Peters, P.S. Robyr, R.W. Bowtell, P. Mansfield, Echo-planar microscopy of porous rocks, Magn. Reson. Imaging 14 (1996) 875–877. J. Hennig, A. Nauerth, H. Friedburg, RARE imaging: a fast imaging method for clinical MR, Magn. Reson. Med. 3 (1986) 823–833. W.P. Rothwell, H.J. Vinegar, Petrophysical applications of NMR imaging, Appl. Opt. 24 (1985) 3969–3972. P. Majors, P. Li, E.J. Peters, NMR imaging of immiscible displacements in porous media, SPE Form. Eval. (1997) 164–169. A.S. Al-Muthana, G.G. Hursan, M. Ma, A. Valori, B. Nicot, P.M. Singer, Wettability as a function of pore size by NMR, in: SCA Conference Paper 31, 2012. S. Emid, J.H.N. Creyghton, High-resolution NMR imaging in solids, Physica B & C 128 (1985) 81–83. B.J. Balcom, R.P. MacGregor, S.D. Beyea, D.P. Green, R.L. Armstrong, T.W. Bremner, Single-point ramped imaging with T1 enhancement (SPRITE), J. Magn. Reson. Ser. A 123 (1996) 131–134. L. Li, H. Han, B.J. Balcom, Spin echo SPI methods for quantitative analysis of fluids in porous media, J. Magn. Reson. 198 (2009) 252–260. Q. Chen, B.J. Balcom, Measurement of rock-core capillary pressure curves using a single-speed centrifuge and one-dimensional magnetic-resonance imaging, J. Chem. Phys. 122 (2005) 214720. Q. Chen, B.J. Balcom, Calillary pressure curve measurement using a single-moderate-speed centrifuge and quantitative magnetic resonance imaging, SCA Conference Paper Number 44, 2005. D.P. Green, J.R. Dick, J. Gardner, B.J. Balcom, B. Zhou, Comparison study of capillary pressure curves obtained using traditional centrifuge and magnetic resonance imaging techniques, SCA Conference Paper Number 30, 2007. D.P. Green, J.R. Dick, M. McAloon, P.F.d.J. Cano-Barrita, J. Burger, B.J. Balcom, Oil/water imbibition and drainage capillary pressure determined by MRI on a wide sampling of rocks, SCA Conference Paper Number 01, 2008. D.P. Green, Extracting pore throat size and relative permeability from MRI based capillary pressure curves, SCA Conference Paper Number 46, 2009. B.A. Baldwin, W.S. Yamanashi, Capillary-pressure determinations from NMR images of centrifuged core plugs: Berea sandstone, Log Anal. (September–October) (1991) 550–556. J.V. Nørgaard, D. Olsen, N. Springer, J. Reffstrup, Capillary pressure curves for low permeability chalk obtained by NMR imaging of core saturation profiles, SPE 30605 (1995) 807–816. C.M. Nielsen, J.K. Larsen, N. Bech, J. Reffstrup, D. Olsen, Determination of saturation functions of tight core samples based on measured saturation profiles, SCA Conference Paper Number 21, 1997.
J. Mitchell et al. / Physics Reports (
)
–
55
[41] E.A. Spinler, B.A. Baldwin, Capillary pressure scanning curves by direct measurement of saturation, SCA Conference Paper Number 05, 1997. [42] K. Romanenko, B.J. Balcom, Permeability mapping in porous media by magnetization prepared centric-scan SPRITE, Expt. Fluid. 50 (2011) 301–312. [43] K. Romanenko, B.J. Balcom, Permeability mapping in naturally heterogeneous sandstone cores by magnetization prepared centric-scan SPRITE, AIChE J. 58 (2012) 3916–3926. http://dx.doi.org/10.1002/aic.13778. [44] M. Bencsik, C. Ramanathan, Direct measurement of porous media local hydrodynamical permeability using gas MRI, Magn. Reson. Imaging 19 (2001) 379–383. [45] P. Coussot, F. Bertrand, B. Herzhaft, Rheological behavior of drilling muds, characterization using MRI visualization, Oil & Gas Sci. Tech. 59 (2009) 23–29. [46] P.T. Callaghan, Rheo-NMR and velocity imaging, Curr. Op. Colloid Interf. Sci. 11 (2006) 13–18. [47] F.R. Wassmuth, K. Green, L. Randall, Details of in situ foam propagation exposed with magnetic resonance imaging, SPE, 59285, 2000, SPE/DOE Improved Oil Recovery Symposium, Tulsa, Oklahoma, 35 April. [48] F.R. Wassmuth, K. Green, L. Randall, Details of in-situ foam propagation exposed with magnetic resonance imaging, SPE Reservoir Eval. & Eng. (April) (2001) 135–145. [49] A. Brautaset, G. Ersland, A. Graue, Determining wettability from in situ pressure and saturation measurements, SCA Conference Paper Number 44, 2010. [50] J. Mitchell, P. Blumler, P.J. McDonald, Spatially resolved nuclear magnetic resonance studies of planar samples, Prog. Nucl. Magn. Reson. Spect. 48 (2006) 161–181. [51] J.M.S. Hutchison, W.A. Edelstein, G. Johnson, A whole-body NMR imaging machine, J. Phys. E: Sci. Instrum. 13 (1980) 947–955. [52] H. Metz, K. Mäder, Benchtop-NMR and MR — a new analytical tool in drug delivery research, Int. J. Pharm. 364 (2008) 170–175. [53] M. Jones, P.S. Aptaker, J. Cox, B.A. Gardiner, P.J. McDonald, A transportable magnetic resonance imaging system for in situ measurements of living trees: the Tree Hugger, J. Magn. Reson. 218 (2010) 133–140. http://dx.doi.org/10.1016/j.jmr.2012.02.019. [54] P.M. Glover, P.S. Aptaker, J.R. Bowler, E. Ciampi, P.J. McDonald, A novel high-gradient permanent magnet for the profiling of planar films and coatings, J. Magn. Reson. 139 (1999) 90–97. [55] J. Perlo, F. Casanova, B. Blümich, Profiles with microscopic resolution by single-sided NMR, J. Magn. Reson. 176 (2005) 64–70. [56] Q. Chen, M.K. Gingras, B.J. Balcom, A magnetic resonance study of pore filling processes during spontaneous imbibition in Berea sandstone, J. Chem. Phys. 119 (2003) 9609–9616. [57] M.D. Hürlimann, L. Venkataramanan, C. Flaum, The diffusion-spin relaxation time distribution function as an experimental probe to characterize fluid mixtures in porous media, J. Chem. Phys. 117 (2002) 10223–10232. [58] Y.Q. Song, L. Venkataramanan, M.D. Hürlimann, M. Flaum, P. Frulla, C. Straley, T1 –T2 correlation spectra obtained using a fast two-dimensional Laplace inversion, J. Magn. Reson. 154 (2002) 261–268. [59] M. Rauschhuber, G. Hirasaki, Determination of saturation profiles via low-field NMR imaging, SCA Conference Paper Number 09, 2009. [60] P. Parasoglou, D.M. Malioutov, A.J. Sederman, J. Rasburn, H. Powell, L.F. Gladden, A. Blake, M.L. Johns, Quantitative single point imaging with compressed sensing, J. Magn. Reson. 201 (2009) 72–80. [61] D.J. Holland, D.M. Malioutov, A. Blake, A.J. Sederman, L.F. Gladden, Reducing data acquisition times in phase-encoded velocity imaging using compressed sensing, J. Magn. Reson. 203 (2010) 236–246. [62] D.J. Holland, J. Mitchell, A. Blake, L.F. Gladden, Particle sizing in porous media using Bayesian magnetic resonance, Phys. Rev. Lett. 110 (2013), article 018001. [63] Y. Chen, B. MacMillan, R.P. MacGregor, B.J. Balcom, Direct detection of hydrocarbon displacement in a model porous soil with magnetic resonance imaging, Anal. Chem. 77 (2005) 1824–1830. [64] D.M. Doddrell, D.T. Pegg, M.R. Bendall, Distortionless enhancement of NMR signals by polarization transfer, J. Magn. Reson. 48 (1982) 323–327. [65] P.N. Tutunjian, H.J. Vinegar, J.A. Ferris, Nuclear magnetic resonance imaging of sodium-23 in cores, SCA Conference Paper Number 11, 1991. [66] P.N. Tutunjian, H.J. Vinegar, J.A. Ferris, Nuclear magnetic resonance imaging of sodium-23 in cores, The Log Anal. (May–June) (1993) 11–19. [67] R. Vold, J. Waugh, M. Klein, D. Phelps, Measurement of spin relaxation in complex systems, J. Chem. Phys. 48 (1968) 3831–3832. [68] T.C. Chandrasekera, J. Mitchell, E.J. Fordham, L.F. Gladden, M.L. Johns, Rapid encoding of T1 with spectral resolution in n-dimensional relaxation correlations, J. Magn. Reson. 194 (2008) 156–161. [69] E.L. Hahn, Nuclear induction due to free larmor precession, Phys. Rev. 77 (1950) 297–298. [70] E.L. Hahn, Spin echoes, Phys. Rev. 80 (1950) 580–594. [71] H. Carr, E. Purcell, Effects of diffusion on free precession in NMR experiments, Phys. Rev. 94 (1954) 630–638. [72] S. Meiboom, D. Gill, Modified spin-echo method for measuring nuclear relaxation times, Rev. Sci. Instrum. 29 (1958) 668–691. [73] R.L. Vold, R.R. Vold, H.E. Simon, Errors in measurements of transverse relaxation rates, J. Magn. Reson. 11 (1973) 283–298. [74] M.D. Hürlimann, D.D. Griffin, Spin dynamics of Carr–Purcell–Meiboom–Gill-like sequences in grossly inhomogeneous B0 and B1 fields and application to NMR well logging, J. Magn. Reson. 143 (2000) 120–135. [75] M.H. Levitt, Spin Dynamics: Basics of Nuclear Magnetic Resonance, John Wiley & Sons Ltd., Chichester, 2001. [76] J. Keeler, Understanding NMR Spectroscopy, John Wiley & Sons Ltd., Chichester, 2005. [77] G. Liu, Y. Li, J. Jonas, Confined geometry effects on reorientational dynamics of molecular liquids in porous silica glasses, J. Chem. Phys. 95 (1991) 6892–6901. [78] S. Godefroy, J.P. Korb, M. Fleury, R.G. Bryant, Surface nuclear magnetic relaxation and dynamics of water and oil in macroporous media, Phys. Rev. E 64 (2001) 021605. [79] K.R. Brownstein, C.E. Tarr, Importance of classical diffusion in NMR studies of water in biological cells, Phys. Rev. A 19 (1979) 2446–2453. [80] S. Davies, K.J. Packer, Pore-size distributions from nuclear magnetic resonance spin-lattice relaxation measurements of fluid-saturated porous solids. I. Theory and simulation, J. Appl. Phys. 67 (1990) 3163–3170. [81] S. Davies, M.Z. Kalam, K.J. Packer, F.O. Zelaya, Pore-size distributions from nuclear magnetic resonance spin-lattice relaxation measurements of fluid saturated porous solids. II. Applications to reservoir core samples, J. Appl. Phys. 67 (1990) 3171–3176. [82] D.F. Allen, A. Boyd, J. Massey, E.J. Fordham, M.O. Amabeoku, W.E. Kenyon, W.B. Ward, The practical application of NMR logging in carbonates: 3 case studies, in: SPLWA 42nd Annual logging symposium, June 17–20, 2001. [83] J.H. Strange, M. Rahman, E.G. Smith, Characterization of porous solids by NMR, Phys. Rev. Lett. 71 (1993) 3589–3591. [84] J.H. Strange, L. Betteridge, M.J.D. Mallett, NATO ASI Series II: Mathematics, Physics and Chemistry. Vol: Magnetic Resonance in Colloid and Interface Science, Kluwer Academic Publishers, Dordrecht, 2002. [85] J. Mitchell, J.B.W. Webber, J.H. Strange, Nuclear magnetic resonance cryoporometry, Phys. Rep. 461 (2008) 1–36. [86] J.E. Tanner, E.D. Stejskal, Restricted self-diffusion of protons in colloidal systems by the pulsed gradient, spin-echo method, J. Chem. Phys. 49 (1968) 1768–1777. [87] K.J. Packer, C. Rees, Pulsed NMR studies of restricted diffusion. I. Droplet size distributions in emulsions, J. Colloid Interf. Sci. 40 (1971) 206–218. [88] P.N. Sen, Time-dependent diffusion coefficient as a probe of geometry, Concept. Magn. Reson. 23A (2004) 1–21. [89] L.L. Latour, R.L. Kleinberg, P.P. Mitra, C.H. Sotak, Pore-size distributions and tortuosity in heterogeneous porous-media, J. Magn. Reson. Ser. A 112 (1995) 83–91. [90] J.D. Wilson, Statistical approach to the solution of 1st kind integral-equations arising in the study of materials and their properties, J. Mater. Sci. 27 (1992) 3911–3924. [91] A.N. Tikhonov, V.Y. Arsenin, Solutions of Ill-Posed Problems, V. H. Winston & Sons, 1977. [92] J. Mitchell, T.C. Chandrasekera, L.F. Gladden, Numerical estimation of relaxation and diffusion distributions in two dimensions, Prog. Nucl. Magn. Reson. Spect. 62 (2012) 34–50.
56
J. Mitchell et al. / Physics Reports (
)
–
[93] G. Wahba, Y.H. Wang, When is the optimal regularization parameter insensitive to the choice of the loss function? Commun. Stat. Theory 19 (1990) 1685–1700. [94] E.D. Stejskal, J.E. Tanner, Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient, J. Chem. Phys. 42 (1965) 288–292. [95] J.E. Tanner, Use of the stimulated echo in NMR diffusion studies, J. Chem. Phys. 52 (1970) 2523–2526. [96] R.L. Kleinberg, W.E. Kenyon, P.P. Mitra, Mechanism of NMR relaxation of fluids in rock, J. Magn. Reson. Ser. A 108 (1994) 206–214. [97] R.M. Cotts, M.J.R. Hoch, T. Sun, J.T. Markert, Pulsed field gradient stimulated echo methods for improved NMR diffusion measurements in heterogeneous systems, J. Magn. Reson. 83 (1989) 252–266. [98] W.S. Price, Pulsed-field gradient nuclear magnetic resonance as a tool for studying translational diffusion: Part I.Basic theory, Concept. Magn. Reson. 9A (1997) 299–336. [99] W.S. Price, Pulsed-field gradient nuclear magnetic resonance as a tool for studying translational diffusion: Part II.Experimental aspects, Concept. Magn. Reson. 10 (1998) 197–237. [100] L.L. Latour, L. Li, C.H. Sotak, Improved PFG stimulated echo method for the measurement of diffusion in inhomogeneous fields., J. Magn. Reson. B 101 (1993) 72–77. [101] M.D. Hürlimann, L. Venkataramanan, Quantitative measurement of two-dimensional distribution functions of diffusion and relaxation in grossly inhomogeneous fields, J. Magn. Reson. 157 (2002) 31–42. [102] M.D. Hürlimann, M. Flaum, L. Venkataramanan, C. Flaum, R. Freedman, G.J. Hirasaki, Diffusion-relaxation distribution functions of sedimentary rocks in different saturation states, Magn. Reson. Imaging 21 (2003) 305–310. [103] J. Mitchell, E.J. Fordham, Emulation of petroleum well-logging D–T2 correlations on a standard benchtop spectrometer, J. Magn. Reson. 212 (2011) 394–401. [104] J. Kärger, W. Heink, The propagator representation of molecular transport in microporous crystallites, J. Magn. Reson. 51 (1983) 1–7. [105] K.J. Packer, J.J. Tessier, The characterization of fluid transport in a porous solid by pulsed gradient stimulated echo NMR, Mol. Phys. 87 (1996) 267–272. [106] J.J. Tessier, K.J. Packer, J.F. Thovert, P.M. Adler, NMR measurements and numerical simulation of fluid transport in porous solids, AIChE J. 43 (1997) 1653–1661. [107] J.J. Tessier, K.J. Packer, The characterization of multiphase fluid transport in a porous solid by pulsed gradient stimulated echo nuclear magnetic resonance, Phys. Fluids 10 (1998) 75–85. [108] K.J. Packer, S. Stapf, J.J. Tessier, R.A. Damion, The characterisation of fluid transport in porous solids by means of pulsed magnetic field gradient NMR, Magn. Reson. Imaging 16 (1998) 463–469. [109] S. Stapf, K.J. Packer, S. Békri, P.M. Adler, Two-dimensional nuclear magnetic resonance measurements and numerical simulations of fluid transport in porous rocks, Phys. Fluids 12 (2000) 566–580. [110] M.L. Johns, A.J. Sederman, L.F. Gladden, A. Wilson, S. Davies, Using MR techniques to probe permeability reduction in rock cores, AIChE J. 49 (2003) 1076–1084. [111] W.M. Holmes, K.J. Packer, Investigation of two phase flow and phase trapping by secondary imbibition within fontainebleau sandstone, Magn. Reson. Imaging 21 (2003) 389–391. [112] W.M. Holmes, K.J. Packer, Investigation ofphase trapping by secondary imbibition within Fontainebleau sandstone, Chem. Eng. Sci. 59 (2004) 2891–2898. [113] U.M. Scheven, D. Verganelakis, R. Harris, M.L. Johns, L.F. Gladden, Quantitative nuclear magnetic resonance measurements of preasymptotic dispersion in flow through porous media, Phys. Fluids 17 (2005) 117107. [114] U.M. Scheven, J.G. Seland, D.G. Cory, NMR-propagator measurements in porous media in the presence of surface relaxation and internal fields, Magn. Reson. Imaging 23 (2005) 363–365. [115] P.M. Singer, G. Leu, E.J. Fordham, P.N. Sen, Low magnetic fields for flow propagators in permeable rocks, J. Magn. Reson. 183 (2006) 167–177. [116] U.M. Scheven, J.P. Crawshaw, V.J. Anderson, R. Harris, M.L. Johns, L.F. Gladden, A cumulant analysis for non-Gaussian displacement distributions in Newtonian and non-Newtonian flows through porous media, Magn. Reson. Imaging 25 (2007) 513–516. [117] J. Mitchell, A.J. Sederman, E.J. Fordham, M.L. Johns, L.F. Gladden, A rapid measurement of flow propagators in porous rocks, J. Magn. Reson. 191 (2008) 267–272. [118] J. Mitchell, D.A. Graf vonder Schulenburg, D.J. Holland, E.J. Fordham, M.L. Johns, L.F. Gladden, Determining NMR flow propagator moments in porous rocks without the influence of relaxation, J. Magn. Reson. 193 (2008) 218–225. [119] J. Mitchell, M.L. Johns, Rapid measurements of diffusion using PFG: developments and applications of the difftrain pulse sequence, Concept. Magn. Reson. 34A (2009) 1–15. [120] W. Zhao, G. Picard, G. Leu, P.M. Singer, Characterization of single-phase flow through carbonate rocks: quantitative comparison of NMR flow propagator measurements with a realistic pore network model, Transp. Porous Med. 81 (2010) 305–315. [121] R. Hussain, T.R.R. Pintelon, J. Mitchell, M.L. Johns, Using NMR displacement measurements to probe CO2 entrapment in porous media, AIChE J. 57 (2010) 1700–1709. [122] E.W. Washburn, P.T. Callaghan, Tracking pore to pore exchange using relaxation exchange spectroscopy, Phys. Rev. Lett. 97 (2006) 175502. [123] M.D. Hürlimann, K.G. Helmer, T.M. deSwiet, P.N. Sen, C.H. Sotak, Spin echoes in a constant gradient and in the presence of simple restriction, J. Magn. Reson. Ser. A 113 (1995) 260–264. [124] P. LeDoussal, P.N. Sen, Decay of nuclear magnetization by diffusion in a parabolic magnetic-field — an exactly solvable model, Phys. Rev. B 46 (1992) 3465–3485. [125] J. Mitchell, T.C. Chandrasekera, L.F. Gladden, Obtaining true transverse relaxation time distributions in high-field NMR measurements of saturated porous media: removing the influence of internal gradients, J. Chem. Phys. 132 (2010) 244705. [126] T.M. deSwiet, P.N. Sen, Decay of nuclear magnetization by bounded diffusion in a constant field gradient, J. Chem. Phys. 100 (1994) 5597–5604. [127] C. Straley, D. Rossini, H.J. Vinegar, C. Morriss, Core analysis by low field NMR, SCA Conference Paper Number 04, 1994. [128] Y.Q. Song, S. Ryu, P.N. Sen, Determining multiple length scales in rocks, Nature 406 (2000) 178–181. [129] Y.Q. Song, Pore sizes and pore connectivity in rocks using the effect of internal field, Magn. Reson. Imaging 19 (2001) 417–421. [130] Y.Q. Song, Using internal magnetic fields to obtain pore size distributions of porous media, Concept. Magn. Reson. 18A (2003) 97–110. [131] M.H.G. Amin, L.D. Hall, R.J. Chorley, T.A. Carpenter, K.S. Richards, B.W. Bache, Magnetic resonance imaging of soil-water phenomena, Magn. Reson. Imaging 12 (1994) 319–321. [132] L.D. Hall, M.H. GaoAmin, E. Dougherty, M. Sanda, J. Votrubova, K.S. Richards, R.J. Chorley, M. Cislerova, MR properties of water in saturated soils and resulting loss of MRI signal in water content detection at 2 tesla, Geoderma 80 (1997) 431–448. [133] A.J. Lucas, M. Peyron, G.K. Pierens, T.A. Carpenter, L.D. Hall, R.C. Stewart, G.F. Potter, D.W. Phelps, A rigorous evaluation of the experimental requirements for performing quantitative porosity measurements by NMR, SPE (July) (1993) 27420. [134] J. Mitchell, J. Staniland, R. Chassagne, E.J. Fordham, Quantitative in-situ enhanced oil recovery monitoring using magnetic resonance, Trans. Porous Med. 94 (2012) 683–706. [135] L. Randall, K. Green, D. Fisher, J. Goldman, T. Prichard, Strategy to reduce uncertainties in core analysis programs through the use of magnetic resonance imaging: a new perspective, SCA Conference Paper Number 36, 2002. [136] L. Venkataramanan, Y.Q. Song, M.D. Hürlimann, Solving Fredholm integrals of the first kind with tensor product structure in 2 and 2.5 dimensions, IEEE Trans. Signal Process. 50 (2002) 1017–1026. [137] N.J. Heaton, Multi-measurement NMR analysis based on maximum entropy, 2005. US patent 6,960,913 B2. [138] J.P. Butler, J.A. Reeds, S.V. Dawson, Estimating solutions of 1st kind integral-equations with nonnegative constraints and optimal smoothing, SIAM J. Numer. Anal. 18 (1981) 381–397.
J. Mitchell et al. / Physics Reports (
)
–
57
[139] P.C. Hansen, Analysis of discrete ill-posed problems by means of the L-curve, SIAM Rev. 34 (1992) 561–580. [140] W. Schoenfelder, H.R. Gläser, I. Mitreiter, F. Stallmach, Two-dimensional NMR relaxometry study of pore space characteristics of carbonate rocks from a Permian aquifer, J. Appl. Geophys. 65 (2008) 21–29. [141] Y.Q. Song, A 2D NMR method to characterize granular structure of dairy products, Prog. Nucl. Magn. Reson. Spect. 55 (2009) 324–334. [142] P.J. McDonald, J.P. Korb, J. Mitchell, L. Monteilhet, Surface relaxation and chemical exchange in hydrating cement pastes: A two-dimensional NMR relaxation study, Phys. Rev. E 72 (2005) 011409. [143] J. Mitchell, M.D. Hürlimann, E.J. Fordham, A rapid measurement of T1 /T2 : the DECPMG sequence, J. Magn. Reson. 200 (2009) 198–206. [144] D. Weber, J. Mitchell, J. McGregor, L.F. Gladden, Comparing strengths of surface interactions for reactants and solvents in porous catalysts using two-dimensional NMR relaxation correlations, J. Phys. Chem. C 113 (2009) 6610–6615. [145] J.G. Seland, M. Bruvold, H. Brurok, P. Jynge, J. Krane, Analyzing equilibrium water exchange between myocardial tissue compartments using dynamical two-dimensional correlation experiments combined with manganese-enhanced relaxography, Magn. Reson. Med. 58 (2007) 674–686. [146] K.M. Song, J. Mitchell, H. Jaffel, L.F. Gladden, Simultaneous monitoring of hydration kinetics, microstructural evolution, and surface interactions in hydrating gypsum plaster in the presence of additives, J. Mater. Sci. 45 (2010) 5282–5290. [147] L. Monteilhet, J.P. Korb, J. Mitchell, P.J. McDonald, Observation of exchange of micropore water in cement pastes by two-dimensional T2 -T2 nuclear magnetic resonance relaxometry, Phys. Rev. E 74 (2006) 061404. [148] J. Mitchell, J.D. Griffith, J.H.P. Collins, A.J. Sederman, L.F. Gladden, M.L. Johns, Validation of NMR relaxation exchange time measurements in porous media, J. Chem. Phys. 127 (2007) 234701. [149] M. Fleury, J. Soualem, Quantitative analysis of diffusional pore coupling from T2 -store-T2 NMR experiments, J. Colloid Interf. Sci. 336 (2009) 250–259. [150] P. Callaghan, I. Furó, Diffusion-diffusion correlation and exchange as a signature for local order and dynamics, J. Chem. Phys. 120 (2004) 4032–4038. [151] S. Arora, D. Horstmann, P. Cherukupalli, J. Edwards, R. Ramamoorthy, T. McDonald, D. Bradley, C. Ayan, J. Zaggas, K. Cig, Single-well in-situ measurement of residual oil saturation after an EOR chemical flood, SPE, 129069, 2010, SPE EOR conference, Muscat, Oman, 11–13 April 2010. [152] G. Leu, E.J. Fordham, M.D. Hürlimann, P. Frulla, Fixed and pulsed gradient diffusion methods in low-field core analysis, Magn. Reson. Imaging 23 (2005) 305–309. [153] M.M. Britton, G.G. Graham, K.J. Packer, Relationships between flow and NMR relaxation of fluids in porous solids, Magn. Reson. Imaging 19 (2001) 325–331. [154] E.J. Fordham, S.J. Gibbs, L.D. Hall, Partially restricted diffusion in a permeable sandstone: observations by stimulated echo PFG NMR, Magn. Reson. Imaging 12 (1994) 279–284. [155] S.J. Gibbs, C.S. Johnson, A PFG NMR experiment for accurate diffusion and flow studies in the presence of eddy currents, J. Magn. Reson. 93 (1991) 395–402. [156] E.M. Haacke, R.W. Brown, M.R. Thompson, R. Venkatesan, Magnetic Resonance Imaging ∼ Physical Principles and Sequence Design, John Wiley & Sons, Inc., 1999. [157] P.G. Morris, Nuclear Magnetic Resonance Imaging in Medicine and Biology, Clarendon Press, Oxford, 1986. [158] A.A. Samo˘ilenko, D.Y. Artemov, L.A. Sibel’dina, Formation of sensitive layer in experiments on NMR subsurface imaging of solids, JEPT Lett. 47 (1988) 417–419. [159] P.J. McDonald, T. Pritchard, S.P. Roberts, Diffusion of water at low saturation levels into sandstone rock plugs measured by broad line magnetic resonance profiling, J. Colloid Interf. Sci. 177 (1996) 439–445. [160] A.J. Lucas, G.K. Pierens, M. Peyron, T.A. Carpenter, L.D. Hall, R.C. Stewart, D. Phelps, G.F. Potter, Quantitative porosity mapping of reservoir rock cores by physically slice selected NMR, SCA Conference Paper Number 01EURO, 1992. [161] G.K. Pierens, M. Peyron, A.J. Lucas, T.A. Carpenter, L.D. Hall, G.F. Potter, R.C. Stewart, D.W. Phelps, Quantitative longitudinal fluid saturation profiles with a slice-selected CPMG sequence, Magn. Reson. Imaging 12 (1994) 323–324. [162] G. Bodenhausen, R. Freedman, G.A. Morris, A simple pulse sequence for selective excitation in fourier transform NMR, J. Magn. Reson. 23 (1976) 171–175. [163] G.A. Morris, R. Freeman, Selective excitation in fourier-transform nuclear magnetic-resonance, J. Magn. Reson. 29 (1978) 433–462. [164] D. Boudot, D. Canet, J. Brondeau, J.C. Boudel, DANTE-Z — a new approach for accurate frequency-selectivity using hard pulses, J. Magn. Reson. 83 (1989) 428–439. [165] C. Roumestand, D. Canet, N. Mahieu, F. Toma, DANTE-Z, an alternative to low power soft pulses — improvement of the selection scheme and applications to multidimensional NMR-studies of proteins, J. Magn. Reson. Ser. A 106 (1994) 168–181. [166] O.V. Petrov, B.J. Balcom, Local T2 distribution measurements with DANTE-Z slice selection, J. Magn. Reson. 215 (2012) 109–114. [167] R. Damadian, Tumor detection by nuclear magnetic resonance, Science 171 (1971) 1151–1153. [168] P. Mansfield, P.K. Grannell, NMR ‘diffraction’ in solids? J. Phys. C: Solid State Phys. 6 (1973) L422–L426. [169] P.C. Lauterbur, Image formation by induced local interactions: examples employing nuclear magnetic resonance, Nature 242 (1973) 190–191. [170] W.P. Rothwell, Nuclear magnetic resonance imaging, Appl. Opt. 24 (1985) 3958–3968. [171] H.J. Vinegar, X-ray CT and NMR imaging of rocks, J. Petrol. Tech. (March) (1986) 257–259. [172] D.A. Doughty, L. Tomutsa, NMR microscopy for fluid imaging at pore scale in reservoir rock, SCA Conference Paper Number 22, 1992. [173] D.A. Doughty, L. Tomutsa, Multinuclear NMR microscopy of two-phase fluid systems in porous rock, Magn. Reson. Imaging 14 (1996) 869–873. [174] W.A. Edelstein, J.M.S. Hutchison, G. Johnson, T. Redpath, Spin warp NMR imaging and applications to human whole-body imaging, Phys. Med. Biol. 25 (1980) 751–756. [175] G. Johnson, J.M.S. Hutchison, T. Redpath, L.M. Eastwood, Improvements in performance time for simultaneous 3-dimensional NMR imaging, J. Magn. Reson. 54 (1983) 374–384. [176] M.P. Enwere, J.S. Archer, NMR imaging for core flood testing, SCA Conference Paper Number 18, 1992. [177] W.A. Edelstein, H.J. Vinegar, P.N. Tutunjian, P.B. Roemer, O.M. Mueller, NMR imaging for core analysis, SPE 18272 (1988) 101–112. 63rd ATC, Houston, Texas, October 2–5, 1988. [178] D.A. Doughty, N.L. Maerefat, Preliminary transformation of an NMR spectrometer into an NMR imager for evaluating fluid content and rock properties of core samples, Log Anal. (March–April) (1989) 78–84. [179] C. Chardaire, J.C. Roussel, NMR imaging of fluid saturation distributions in core samples using a high magnetic field, SCA Conference Paper Number 15EURO, 1990. [180] P.A. Osment, K.J. Packer, M.J. Taylor, J.J. Attard, T.A. Carpenter, L.D. Hall, N.J. Herrod, S.J. Doran, NMR imaging of fluids in porous solids, Philos. Trans. R. Soc. Lond. A Math. Phys. Eng. Sci. 333 (1990) 441–452. [181] M.A. Robinson, Magnetic resonance imaging for the measurement of oil core porosity distributions, SPE (July) (1991) 23582. [182] J.W. Gleeson, D.E. Woessner, Three-dimensional and flow-weighted NMR imaging of pore connectivity in a limestone, Magn. Reson. Imaging 9 (1991) 879–884. [183] M.R. Merrill, Porosity measurements in natural porous rocks using magnetic resonance imaging, Appl. Magn. Reson. 5 (1993) 307–321. [184] M.R. Merrill, Local velocity and porosity measurements inside casper sandstone using MRI, AIChE J. 40 (1994) 1262–1267. [185] K. Mirotchnik, P. Kubika, L. Randall, A. Starosud, A. Allsopp, K. amdKantzas, Determination of mud invasion characteristics of sandstone reservoirs using a combination of advanced core analysis techniques, SCA Conference Paper Number 15, 1998. [186] P. Mansfield, Multi-planar image formation using NMR spin echoes, J. Phys. C: Solid State Phys. 10 (1977) L55–L58. [187] P. Mansfield, R. Bowtell, S. Blackband, D.N. Guilfoyle, Magnetic resonance imaging: applications of novel methods in studies of porous media, Magn. Reson. Imaging 10 (1992) 741–746. [188] M.D. Mantle, A.J. Sederman, Dynamic mri in chemical process and reaction engineering, Prog. Nucl. Magn. Reson. Spect. 43 (2003) 3–60. [189] D.N. Guilfoyle, P. Mansfield, Fluid flow measurement in porous media by echo-planar imaging, J. Magn. Reson. 97 (1992) 342–358.
58
J. Mitchell et al. / Physics Reports (
)
–
[190] G. Maddinelli, C. Carati, A. Brancolini, E. Rossi, E. Causin, NMR imaging of gelation processes inside rock cores, SPE 25219 (1993) International Symposium on Oilfield Chemistry, New Orleans, LA, USA, March 2–5. [191] E.H. Isaaks, R.M. Srivastava, An Introduction to Applied Geostatistics, Oxford University Press, Oxford, 1989. [192] S. Bansal, M.R. Merrill, H.A. Deans, Application of the semivariogram to MRI data: a useful approach for the characterization of natural porous rocks, Appl. Magn. Reson. 4 (1993) 241–259. [193] A.E. Pomerantz, P. Tilke, Y.Q. Song, Inverting MRI measurements to heterogeneity spectra, J. Magn. Reson. 193 (2008) 243–250. [194] S. Baraka-Lokmane, B.T. Ngwenya, I.G. Main, Benefit of complementary methods for characterizing sandstone cores, SCA Conference Paper Number 57, 2007. [195] S. Baraka-Lokmane, I.G. Main, B.T. Ngwenya, S.C. Elphick, Application of complementary methods for more robust characterization of sandstone cores, Marine Petrol. Geo. 26 (2009) 39–56. [196] A. Bersani, B.B. Bam, F. Radaelli, E. Rossi, Improved methodology for the characterization of complex vuggy carbonate, SCA Conference Paper Number 10, 2009. [197] G.C. Borgia, V. Bortolotti, P. Fantazzini, Combined magnetic resonance relaxation and imaging for quantitative determination of matrix and vugular porosity in carbonate cores, SPE 49294 (1998) Annual Technical Conference and Exhibition, New Orleans, Louisiana, 27–30 September. [198] V. Bortolotti, P. Macini, E. Mesini, F. Srisuriyachai, P. Fantazzini, M. Gombia, Combined spatially resolved and non-resolved 1 H-NMR relaxation analysis to assess and monitor wettability reversal in carbonate rocks, IPTC 13443 (2009) IPTC, Doha, Qatar, 7–9 December, 2009. [199] V. Bortolotti, P. Macini, E. Mesini, F. Srisuriyachai, P. Fantazzini, M. Gombia, Probing wettability reversal in carbonatic rocks by spatially-resolved and non-resolved 1H-NMR relaxation analysis, SPE 133937 (2010) Annual Technical Conference and Exhibition, Florence, Italy, 19–22 September. [200] G. Maddinelli, A. Camanzi, C. Carati, A. Brancolini, Heterogeneous systems: an integrated approach based on different imaging techniques, SCA Conference Paper Number 06, 1994. [201] C.H. Arns, An analysis of NMR-permeability scaling rules by numerical MRI, in: SPWLA 48th Annual Logging Symposium, June 3–6, 2007. [202] G.S. Padhy, M.A. Ioannidis, C. Lemaire, Special core analysis studies in vuggy porous media of controlled microstructure, SCA Conference Paper Number 76, 2005. [203] G.S. Padhy, C. Lemaire, M.A. Ioannidis, Magnetic resonance methods for the characterization of the pore space in vuggy carbonates, SCA Conference Paper Number 30, 2006. [204] N. Bona, B. Bam, M. Pirrone, E. Rossi, Use of a new impedance cell and 3D MRI to obtain fast and accurate resistivity index measurements from a single centrifuge step, SCA Conference Paper Number 22, 2011. [205] S. Chen, F. Qin, K.H. Kim, A.T. Watson, NMR imaging of multiphase flow in porous media, SPE 24760 (1992) 67th Annual Technical Conference and Exhibition of the Society of Petroleum Engineers held in Washington, DC, October 4–7. [206] S. Chen, F. Qin, K.H. Kim, A.T. Watson, NMR imaging of multiphase flow in porous media, AIChE J. 39 (1993) 925–934. [207] A. Brancolini, A. Cominelli, R. Kulkarni, A.T. Watson, Spatial distribution of petrophysical parameters on a core scale using magnetic resonance imaging, in: SPWLA 38th Annual Logging Symposium, June 15-18, 1997. [208] D. Mattiello, M. Balzarini, L. Ferraccioli, A. Brancolini, Calculation of constituent porosity in a dual-porosity matrix: MRI and image analysis integration, SCA Conference Paper Number 06, 1997. [209] M.A. Robinson, Effects of RF attenuation in magnetic resonance imaging on quantitative measurements in oil cores, J. Magn. Reson. 95 (1991) 574–580. [210] M.A. Horsfield, E.J. Fordham, C. Hall, L.D. Hall, 1 H NMR imaging studies of filtration in colloidal suspensions, J. Magn. Reson. 81 (1989) 593–596. [211] J.J. Attard, S.J. Doran, N.J. Herrod, T.A. Carpenter, L.D. Hall, Quantitative NMR spin-lattice-relaxation imaging of brine in sandstone reservoir cores, J. Magn. Reson. 96 (1992) 514–525. [212] Y.Y. Chen, L.P. Hughes, L.F. Gladden, M.D. Mantle, Quantitative ultra-fast MRI of HPMC swelling and dissolution, J. Pharm. Sci. 99 (2010) 3462–3472. [213] Q. Chen, A.E. Marble, B.G. Colpitts, B.J. Balcom, The internal magnetic field distribution, and single exponential magnetic resonance free induction decay, in rocks, J. Magn. Reson. 175 (2005) 300–308. [214] Q. Chen, H. Hickey, B.J. Balcom, A single-shot method for determining drainage and imbibition capillary pressure curves, SCA Conference Paper Number 12, 2006. [215] A. Belonogov, F. Marica, A. Lawfield, K. Butler, Q. Chen, M.K. Gingras, B.J. Balcom, Magnetic resonance imaging of porosity heterogeneity in bioturbated sandstone from the white rose reservoir, atlantic canada, SCA Conference Paper Number 60, 2005. [216] L.B. Romero-Zerón, S. Ongsurakul, L. Li, B.J. Balcom, Visualization of the effect of porous media wettability on polymer flooding performance through unconsolidated porous media using magnetic resonance imaging, Petrol. Sci. Tech. 28 (2010) 52–67. [217] L.B. Romero-Zerón, S. Ongsurakul, L. Li, B.J. Balcom, Magnetic resonance imaging of phase trapping and in situ permeability modification in unconsolidated porous media, Petrol. Sci. Tech. 28 (2010) 262–276. [218] L. Li, F. Marica, Q. Chen, B. MacMillan, B.J. Balcom, Quantitative discrimination of water and hydrocarbons in porous media by magnetisation prepared centric-scan SPRITE, J. Magn. Reson. 186 (2007) 282–292. [219] F.R. Rack, B.J. Balcom, R.P. MacGregor, R.L. Armstrong, Magnetic resonance imaging of the lake agassiz-lake winnipeg transition, J. Paleolimnology 19 (1997) 255–264. [220] M.K. Gringas, B. MacMillan, B.J. Balcom, Visualizing the internal physical characteristics of carbonate sediments with magnetic resonance imaging and petrography, Bull. Can. Petrol. Geo. 50 (2002) 363–369. [221] F. Marica, Q. Chen, A. Hamilton, C. Hall, T. Al, B.J. Balcom, Spatially resolved measurement of rock core porosity, J. Magn. Reson. 178 (2006) 136–141. [222] O.V. Petrov, B.J. Balcom, Two-dimensional T2 distribution mapping in porous solids with phase encode MRI, J. Magn. Reson. 212 (2011) 102–108. [223] J. Kaffanke, T. Dierkes, S. Romanzetti, M. Halse, J. Rioux, M.O. Leach, B.J. Balcom, N.J. Shah, Application of the chirp z-transform to MRI data, J. Magn. Reson. 178 (2006) 121–128. [224] J. Rioux, M. Halse, E. Aubanel, B.J. Balcom, J. Kaffanke, S. Romanzetti, T. Dierkes, N.J. Shah, An accurate nonuniform Fourier transform for SPRITE magnetic resonance imaging data, ACM Trans. Math. Software 33 (2007) 1–21. [225] J.C. García-Naranjo, P.M. Glover, F. Marica, B.J. Balcom, Variable bandwidth filtering for magnetic resonance imaging with pure phase encoding, J. Magn. Reson. 202 (2010) 234–238. [226] S.D. Beyea, B.J. Balcom, I.V. Mastikhin, T.W. Bremner, R.L. Armstrong, P.E. Grattan-Bellew, Imaging of heterogeneous materials with a turbo spin echo single-point imaging technique, J. Magn. Reson. 144 (2000) 255–265. [227] S.D. Beyea, T.W. Bremner, B.J. Balcom, Minimization of diffusive attenuation in T2 -weighted NMR images of porous solids using turboSPI, Solid State Nucl. Magn. Reson. 29 (2006) 267–271. [228] O.V. Petrov, G. Ersland, B.J. Balcom, T2 distribution mapping profiles with phase-encode MRI, J. Magn. Reson. 209 (2011) 39–46. [229] T. Gullion, D.B. Baker, M.S. Conradi, New, compensated Carr–Purcell sequences, J. Magn. Reson. 89 (1990) 479–484. [230] A.A. Maudsley, S.K. Hilal, W.H. Perman, H.E. Simon, Spatially resolved high resolution spectroscopy by four-dimensional NMR, J. Magn. Reson. 51 (1983) 147–152. [231] R. Kimmich, D. Hoepfel, Volume-selective multipulse spin-echo spectroscopy, J. Magn. Reson. 72 (1987) 379–384. [232] J.M. Dereppe, C. Moreaux, K. Schenker, 2D spin-echo and 3D chemical-shift-imaging techniques for analysis of oil-water replacement in limestone, J. Magn. Reson. 91 (1991) 596–603. [233] J.M. Dereppe, C. Moreaux, K. Schenker, Chemical shift imaging of fluid filled porous rocks, Magn. Reson. Imaging 9 (1991) 809–813. [234] J. Dechter, R. Komoroski, S. Ramaprasad, NMR chemical shift selective imaging of individual fluids in sandstone and dolomite cores, SCA Conference Paper Number 03, 1989. [235] J. Dechter, R. Komoroski, S. Ramaprasad, Use of presaturation for chemical-shift-selective imaging of individual fluids in sandstone and carbonate cores, J. Magn. Reson. 93 (1991) 142–150.
J. Mitchell et al. / Physics Reports (
)
–
59
[236] W.T. Dixon, Simple proton spectroscopic imaging, Radiology 153 (1984) 189–194. [237] P.D. Majors, J.L. Smith, F.S. Kovarik, E. Fukushima, NMR spectroscopic imaging of oil displacement in dolomite, J. Magn. Reson. 89 (1990) 470–478. [238] M.A. Horsfield, C. Hall, L.D. Hall, Two-species chemical-shift imaging using prior knowledge and estimation theory. Application to rock cores, J. Magn. Reson. 87 (1990) 319–330. [239] E.J. Fordham, M.A. Horsfield, C. Hall, L.D. Hall, Low-contrast secondary imbibition in long rock cores, Magn. Reson. Imaging 9 (1991) 803–808. [240] C.T. Chang, C.M. Edwards, Proton MR two-component chemical shift imaging of fluids in porous media, SCA Conference Paper Number 14, 1991. [241] C.T. Chang, C.M. Edwards, Proton MR two-component chemical shift imaging of fluids in porous media, The Log Anal. (November–December) (1993) 20–28. [242] D.N. Saraf, I. Fatt, Three phase relative permeability measurement using a nuclear magnetic resonance technique for estimating fluid saturation, SPE J. (September) (1967) 235–242. [243] G. Maddinelli, A. Brancolini, MRI as a tool for the study of waterflooding processes in heterogeneous cores, Magn. Reson. Imaging 14 (1996) 915–917. [244] S. Mahmood, D.A. Doughty, L. Tomutsa, M.M. Honarpour, Pore level fluid imaging using high resolution nuclear magnetic resonance imaging and thin slab micromodels, SCA Conference Paper Number 24, 1990. [245] S.N. Sarkar, J.J. Dechter, R. Komoroski, Multinuclear NMR imaging of fluid phases in Berea sandstone, J. Magn. Reson. Ser. A 102 (1993) 314–317. [246] L.D. Hall, V. Rajanayagam, Thin-slice, chemical-shift imaging of oil and water in sandstone rock at 80 MHz, J. Magn. Reson. 74 (1987) 139–146. [247] J.L.A. Williams, G. Maddinelli, D.G. Taylor, Selective magnetic resonance imaging of chemicals in sandstone cores during flow, J. Magn. Reson. Ser. A 109 (1994) 124–128. [248] S. Blackband, P. Mansfield, J.R. Barnes, A.D.H. Clague, S.A. Rice, Discrimination of crude oil and water in sand and in bore cores with NMR imaging, SPE Form. Eval. (February) (1986) 31–34. [249] G. Guillot, C. Chardaire-Rivière, S. Bobroff, A. LeRoux, J.C. Roussel, L. Cuiec, Characterisation of wetting heterogeneities in sandstone rocks by MRI, Magn. Reson. Imaging 12 (1994) 365–368. [250] S. Davies, A. Hardwick, D. Roberts, K. Spowage, K.J. Packer, Quantification of oil and water in preserved reservoir rock by NMR spectroscopy and imaging, Magn. Reson. Imaging 12 (1994) 349–353. [251] B.A. Baldwin, W.S. Yamanashi, NMR imaging of fluid dynamics in reservoir core, Magn. Reson. Imaging 6 (1988) 493–500. [252] A.E. Fischer, B.J. Balcom, E.J. Fordham, T.A. Carpenter, L.D. Hall, A fast inversion recovery NMR imaging technique for mapping two-dimensional tracer diffusion and dispersion in heterogeneous media, J. Phys. D: Appl. Phys. 28 (1995) 384–397. [253] R.L. Kleinberg, H.J. Vinegar, NMR properties of reservoir fluids, The Log Anal. (November–December) (1996) 20–32. [254] J. Mitchell, J. Staniland, E.J. Fordham, Paramagnetic doping agents in magnetic resonance studies of oil recovery: a fresh perspective, Petrophysics (2013) (submitted for publication). [255] M.A. Mahmoud, H.A. Nasr-El-Din, C.A. DeWolf, A.K. Alex, Effect of lithology on the flow of chelating agents in porous media during matrix acid treatments, SPE, 140149, 2011, Production and Operations Symposium, Oklahoma, USA, 27–29 March. [256] B. Sun, K.J. Dunn, G.A. LaTorraca, D.M. Wilson, NMR imaging with diffusion and relaxation, SCA Conference Paper Number 24, 2003. [257] T.A. Carpenter, E.S. Davies, C. Hall, L.D. Hall, W.D. Hoff, M.A. Wilson, Capillary water migration in rock: process and material properties examined by NMR imaging, Mater. Struct. 26 (1993) 286–292. [258] M. Gombia, V. Bortolotti, R.J.S. Brown, M. Camaiti, P. Fantazzini, Models of water imbibition in untreated and treated porous media validated by quantitative magnetic resonance imaging, J. Appl. Phys. 103 (2008) 094913. [259] E.J. Fordham, L.D. Hall, T.S. Ramakrishnan, M.R. Sharpe, C. Hall, Saturation gradients in drainage of porous media: NMR imaging measurements, AIChE J. 39 (1993) 1431–1443. [260] S. Wickramathilaka, J.J. Howard, J.C. Stevens, N.R. Morrow, G. Mason, Å. Haugen, A. Graue, M.A. Fernø, Magnetic resonance imaging of oil recovery during spontaneous imbibition from cores with two-ends open boundary condition, SCA Conference Paper Number 24, 2011. [261] B.A. Baldwin, E.A. Spinler, In-situ saturation development during spontaneous imbibition, SCA Conference Paper Number 22, 1999. [262] G. Guillot, A. Trokiner, L. Darrasse, H. Saint-Jaimes, Drying of a porous rock monitored by NMR imaging, J. Phys. D: Appl. Phys. 22 (1989) 1646–1649. [263] C. Straley, D. Rossini, L.M. Schwartz, Particle filtration in sandstone cores: application of a chemical shift imaging technique, SCA Conference Paper Number 17, 1992. [264] C. Straley, D. Rossini, L.M. Schwartz, M.E. Stromski, M. Hrovat, S. Patz, Chemical shift imaging of particle filtration in sandstone cores, Magn. Reson. Imaging 12 (1994) 313–315. [265] C. Straley, D. Rossini, L.M. Schwartz, M.E. Stromski, M. Hrovat, S. Patz, Particle filtration in sandstone cores: a novel application of chemical shift magnetic resonance imaging techniques, The Log Anal. (March–April) (1995) 42–51. [266] C.H. vander Zwaag, F. Stallmach, J.E. Hanssen, Application of NMR- and CT-analytical methods to assess the formation damage potential of drilling fluids, SCA Conference Paper Number 13, 1998. [267] C.H. vander Zwaag, R. Crossley, G. van Graas, Potential effects of current and future drilling fluid systems on core analysis, SCA Conference Paper Number 16, 2000. [268] C.H. vander Zwaag, Benchmarking the formation damage of drilling fluids, SPE 86544 (2004) International Symposium on Formation Damage Control, Lafayette, Louisiana, USA, 18–20 March. [269] R.B. Gharbi, N. Smaoui, E.J. Peters, Unstable EOR displacements and their prediction using the Karhunen–Loéve (K–L) decomposition, SPE 39882 (1998) International Petroleum Conference and Exhibition, Vilkhannosa, Mexico, 3–5 March. [270] E. Aspenes, G. Ersland, A. Graue, J. Stevens, B.A. Baldwin, Wetting phase bridges establish capillary continuity across open fractures and increase oil recovery in mixed-wet fractured chalk, Trans. Porous Med. 74 (2008) 35–47. [271] B.A. Baldwin, W.S. Yamanashi, NMR imaging of fluid saturation distributions in cores, SCA Conference Paper Number 14, 1987. [272] B.A. Baldwin, W.S. Yamanashi, Detecting fluid movement and isolation in reservoir core with medical NMR imaging techniques, SPE Reservoir Engineering (May) (1989) 207–212. [273] J.L.A. Williams, M.P. Enwere, J.S. Archer, D.G. Taylor, Applications of magnetic resonance imaging in special core analysis studies, SCA Conference Paper Number 16EURO, 1990. [274] J.L.A. Williams, D.G. Taylor, G. Maddinelli, M.P. Enwere, J.S. Archer, Visualisation of fluid displacement in rock cores by NMR imaging, Magn. Reson. Imaging 9 (1991) 767–773. [275] M.P. Enwere, J.S. Archer, NMR imaging for water/oil displacement in cores under viscous-capillary force control, SPE/DOE 24166 (1992) 99–104. SPE/DOE 8th Symposium on EOR, Tulsa, Oklahomo, April 22–24, 1992. [276] D. Olsen, S. Topp, A. Stensgaard, J.V. Nørgaard, J. Reffstrup, Quantitative 1D saturation profiles on chalk by NMR, Magn. Reson. Imaging 14 (1996) 847–851. [277] K.H. Kim, S. Chen, F. Qin, A.T. Watson, Use of NMR imaging for determining fluid saturation distributions during multiphase displacement in porous media, SCA Conference Paper Number 19, 1992. [278] H. Riskedal, L. Tipura, J. Howard, A. Graue, NMR monitoring of spontaneous brine imbibition in carbonates, SCA Conference Paper 13, 2008. [279] N. Loahardjo, N.R. Morrow, J. Stevens, J. Howard, Nuclear magnetic resonance imaging: application to determination of saturation changes in a sandstone core by sequential waterflooding, SCA Conference Paper Number 16, 2010. [280] R.N. Kulkarni, A.T. Watson, J.E. Nordtvedt, A. Brancolini, O. Johnsen, Estimation of multiphase flow functions from dynamic displacement applications of NMR imaging, SPE 36855 (1996) European Petroleum Conference, MiIan, Italy, 22–24 October. [281] R. Kulkarni, A.T. Watson, J.E. Nordtvedt, Estimation of porous media flow functions using NMR imaging data, Magn. Reson. Imaging 16 (1998) 707–709. [282] J.L.A. Williams, D.G. Taylor, Measurements of viscosity and permeability of two phase miscible fluid flow in rock cores, Magn. Reson. Imaging 12 (1994) 317–318.
60
J. Mitchell et al. / Physics Reports (
)
–
[283] N. Bech, D. Olsen, C.M. Nielsen, Determination of oil/water saturation functions of chalk core plugs from two-phase flow experiments, SPE Reservoir Eval. & Eng. 3 (2000) 50–59. [284] C. Laroche, J. Kamath, F. Nakagawa, Y. Yortsos, Gas relative permeability reduction due to waterblocking — a laboratory based analysis, SCA Conference Paper Number 38, 2000. [285] M.A. Barrufet, J.M. Perez, S.S. Mondava, L. Ali, S.W. Poston, MRI studies of permeability control using biopolymers, SPE 24809 (1992) 67h Annual Technical Conference and Exhibition, Washington, DC, October 4–7. [286] M.A. Barrufet, R.W. Flumerfelt, M.P. Walsh, A.T. Watson, Development of nuclear magnetic resonance imaging/spectroscopy for improved petroleum recovery, U. S. Department of Energy, DOE/BC/14446-10, 2004. [287] C. Chardaire-Rivière, J.C. Roussel, Use of a high magnetic field to visualize fluids in porous media by MRI, SCA Conference Paper Number 12, 1991. [288] S. Nicula, A. Brancolini, A. Cominelli, S. Mantica, Blending magnetic resonance imaging and numerical simulation, SCA Conference Paper Number 19, 1998. [289] E. Zuluaga, P.D. Majors, E.J. Peters, A simulation approach to validate petrophysical data from NMR imaging, SPE J. (March) (2002) 35–39. [290] S. Chen, X. Yao, J. Qiao, A.T. Watson, NMRI characterization of fractures and multiphase transport in fractured porous media, SPE 28369 (1994) 69th Annual Technical Conference and Exhibition, New Orleans, LA, USA, 25–28 September. [291] S. Chen, X. Yao, J. Qiao, A.T. Watson, Characterization of fractures and multiphase flow in fractured permeable rocks using NMR imaging techniques, SCA Conference Paper Number 02, 1994. [292] C.M. Nielsen, D. Olsen, N. Bech, Imbibition processes in fractured chalk core plugs with connate water mobilization, SPE 63226 (2000) Annual Technical Conference and Exhibition, Dallas, Texas, 14 October. [293] M.A. Fernø, Å. Haugen, A. Graue, J.J. Howard, The significance of wettability and fracture properties on oil recovery efficiency in fractured carbonates, SCA Conference Paper Number 22, 2008. [294] D. Olsen, C.M. Nielsen, N. Bech, Variation in fracture permeability in chalk core plugs, SPE 60302 (2000) Rocky Mountain Regional/Low Permeability Reservoirs Symposium and Exhibition, Denver, Colorado, 12–15 March. [295] D.P. Tobola, B.A. Baldwin, H.E. Farrell, Geological description and flow characterization of a north sea chalk containing stylolites and a gouge-filled fracture, SCA Conference Paper Number 32, 1998. [296] A.T.A. Kumar, P. Majors, W. Rossen, Measurement of aperture and multiphase flow in fracture with NMR imaging, SPE Form. Eval. (June) (1997) 101–107. [297] A. Graue, B.A. Baldwin, E. Aspenes, J. Stevens, D.P. Tobola, D.R. Zornes, Complementary imaging techniques applying NTI and MRI determined wettability effects on oil recovery mechanisms in fractured reservoirs, SCA Conference Paper Number 08, 2003. [298] G. Ersland, M.A. Ferno, A. Graue, B.A. Baldwin, J. Stevens, Complementary imaging of oil recovery mechanisms in fractured reservoirs, Chem. Eng. J. 158 (2010) 32–38. [299] S. Baraka-Lokmane, G. Teutsch, I.G. Main, Influence of open and sealed fractures on fluid flow and water saturation in sandstone cores using magnetic resonance imaging, Geophys. J. Int. 147 (2001) 263–271. [300] A. Graue, E. Aspenes, R.W. Moe, B.A. Baldwin, A. Moradi, J. Stevens, D.P. Tobola, MRI tomography of saturation development in fractures during waterfloods at various wettability conditions, SPE 71506 (2001) Annual Technical Conference and Exhibition, New Orleans, Louisiana, 30 September–3 October. [301] E. Aspenes, A. Graue, B.A. Baldwin, A. Moradi, J. Stevens, D.P. Tobola, Fluid flow in fractures visualized by MRI during waterfloods at various wettability conditions — enmphasis on fracture width and flow rate, SPE 77338 (2002) Annual Technical Conference and Exhibition, San Antonio, Texas, 29 September–2 October. [302] G. Ersland, A. Graue, R.W. Moe, E. Aspenes, T. Skar, S. Berg, Impact of deformation bands on fluid flow and oil recovery, in: SCA Conference Paper Number 45, 2010. [303] J. Stevens, B.A. Baldwin, A. Graue, G. Ersland, J. Husebø, J.J. Howard, Measurement of hydrate formation in sandstone, Petrophysics 49 (2008) 67–73. [304] R.D. Hazlett, J.W. Gleeson, H. Laali, R. Navarro, NMR imaging application to carbon dioxide miscible flooding in west texas carbonates, SCA Conference Paper Number 11, 1993. [305] T. Suekane, N. Furukawa, S. Tsushima, S. Hirai, M. Kiyota, Application of MRI in the measurement of two-phase flow of supercritical CO2 and water in porous rocks, J. Porous Med. 12 (2009) 143–154. [306] A. Brautaset, G. Ersland, A. Graue, J. Stevens, J. Howard, Using MRI to study in situ oil recovery during CO2 injection in carbonates, SCA Conference Paper Number 41, 2008. [307] Y. Wang, T. Luce, C. Ishizawa, M. Shuck, K. Smith, H. Ott, M. Appel, Halite precipitation and permeability assessment during supercritical CO2 core flood, SCA Conference Paper Number 18, 2010. [308] D.P. Tobola, B.A. Baldwin, Capillary continuity in fractured chalk systems: an experimental study, SCA Conference Paper Number 26, 1995. [309] L. Moudrakovski, G. McLaurin, C. Ratcliffe, J. Ripmeester, Methane and carbon dioxide hydrate formation in water droplets: spatially resolved measurements from magnetic resonance microimaging, J. Chem. Phys. B 108 (2004) 17591–17595. [310] A. Graue, E. Aspenes, B. Kvamme, B.A. Baldwin, J. Stevens, D.P. Tobola, D.R. Zornes, Monitoring the formation and dissociation of gas hydrate in reservoir rock using magnetic resonance imaging, SCA Conference Paper Number 55, 2003. [311] J. Stevens, B.A. Baldwin, A. Graue, G. Ersland, J. Husebø, J.J. Howard, Measurements of hydrate formation in sandstone, SCA Conference Paper Number 36, 2007. [312] A. Graue, B. Kvamme, B.A. Baldwin, J. Stevens, J. Howard, E. Aspenes, G. Ersland, J. Husebø, D. Zomes, MRI visualization of spontaneous methane production from hydrates in sandstone core plugs when exposed to CO2 , SPE J. (June) (2008) 146–152. [313] P. Mansfield, B. Issa, Fluid transport in porous rocks. I. EPI studies and a stochastic model of flow, J. Magn. Reson. Ser. A 122 (1996) 137–148. [314] S. Davies, A. Hardwick, K. Spowage, K.J. Packer, Fluid velocity imaging of reservoir core samples, Magn. Reson. Imaging 12 (1994) 265–268. [315] P. Mansfield, B. Issa, Studies of fluid transport in porous rocks by echo-planar MRI, Magn. Reson. Imaging 12 (1994) 275–278. [316] P. Mansfield, B. Issa, Fluid transport in porous rocks. II. Hydrodynamic model of flow and intervoxel coupling, J. Magn. Reson. Ser. A 122 (1996) 149–156. [317] M.A. Al-Mugheiry, B. Issa, P. Mansfield, Imaging fluid movements through sandstones, sands and model glass-bead packs using fast NMR imaging techniques, SPE 68106 (2001) Middle East Oil Show, Bahrain, 17–20 March. [318] R.A. Waggoner, E. Fukushima, Velocity distribution of slow fluid flows in bentheimer sandstone: an NMRI and propagator study, Magn. Reson. Imaging 14 (1996) 1085–1091. [319] C.M. Edwards, C.T. Chang, S.N. Sarkar, The measurement of fluid velocities in porous media, SCA Conference Paper Number 10, 1993. [320] M.R. Merrill, Z. Jin, Velocity measurements in natural porous rocks, Magn. Reson. Imaging 12 (1994) 345–348. [321] G.L. Hassler, E. Brunner, Measurement of capillary pressures in small samples, Petrol. Trans. AIME 160 (1945) 114–123. [322] P. Forbes, Simple and accurate methods for converting centrifuge data into drainage and imbibition capillary pressure curves, Log Anal. 35 (1994) 31–52. [323] D.W. Ruth, Z.A. Chen, Measurement and interpretation of centrifuge capillary pressure curves — the SCA survey data, Log Anal. 36 (1995) 21–32. [324] E.W. Washburn, Note on a method of determining the distribution of pore sizes in a porous material, Proc. Natl. Acad. Sci. 7 (1921) 115–116. [325] B.A. Baldwin, W.S. Yamanashi, NMR imaging of fluid saturation distributions in cores, Log Anal. (September–October) (1991) 536–549. [326] M.C. Leverett, Capillary behavior in porous solids, Trans. AIME 142 (1941) 152–169. Tulsa Meeting, October 1940. [327] Z.A. Chen, D.W. Ruth, The effect of gravity degredation on low-speed centrifuge capillary pressure data, AIChE J. 41 (1995) 469–480. [328] B.A. Baldwin, E.A. Spinler, A direct method for simultaneously determining positive and negative capillary pressure curves in reservoir rock, J. Petrol. Sci. Eng. 20 (1998) 161–165.
J. Mitchell et al. / Physics Reports (
)
–
61
[329] E.A. Spinler, B.A. Baldwin, A. Graue, Simultaneous measurement of multiple capillary pressure curves from wettability and rock property variations within single rock plugs, SCA Conference Paper Number 57, 1999. [330] T. Bognø, E. Aspenes, A. Graue, E.A. Spinler, D.P. Tobola, Comparison of capillary pressure measurements at various wettabilities using the direct measurement of saturation method and conventional centrifuge techniques, SCA Conference Paper Number 28, 2001. [331] B.A. Baldwin, D. Green, J.C. Stevens, J.J. Howard, Sensitivity study of variables associated with magnetic resonance imaging capillary pressure measurements in porous rocks, in: SCA Conference Paper Number 32, 2009. [332] E.A. Spinler, D.R. Maloney, Application of linear X-ray analysis using absorption coefficients for direct determination of in situ core saturation for pc measurement, SCA Conference Paper Number 31, 2001. [333] E.C. Donaldson, R.D. Thomas, P.B. Lorenz, Wettability determination and its effect on recovery efficiency, SPE J. (March) (1969) 13–20. [334] J.V. Nørgaard, D. Olsen, J. Reffstrup, N. Springer, Capillary-pressure curves for low-permeability chalk obtained by nuclear magnetic resonance imaging of core-saturation profiles, SPE Reservoir Eval. & Eng. 2 (1999) 141–148. [335] C. Lindsay, C. Cornwell, D. Green, Acquisition of core capillary pressure data by in-situ saturation monitoring — a comparative evaluation, SCA Conference Paper Number 20, 2009. [336] H. Ritter, L. Drake, Pore-size distribution in porous materials, Industrial Eng. Chem. 17 (1945) 782–786. [337] J. Mitchell, J. Staniland, A. Wilson, A. Howe, A. Clarke, E.J. Fordham, J. Edwards, R. Faber, R. Bouwmeester, Magnetic resonance imaging of chemical EOR in core to complement field pilot studies, SCA Conference Paper Number 34, 2012. [338] J. Mitchell, J. Staniland, R. Chassagne, K. Mogensen, S. Frank, E.J. Fordham, Continuous in-situ monitoring of secondary oil recovery from a microporous limestone using low field magnetic resonance relaxation mapping, J. Petrol. Sci. Eng. (2013) (submitted for publication). [339] J. Mitchell, J. Edwards, E.J. Fordham, J. Staniland, R. Chassagne, P. Cherukupalli, O. Wilson, R. Faber, R. Bouwmeester, Quantitative remaining oil interpretation using magnetic resonance: from the laboratory to the pilot, SPE 154704 (2012) EOR Conference at Oil and Gas West Asia held in Muscat, Oman, 16–18 April 2012. [340] D. Green, D. Veselinovic, B.J. Balcom, F. Marica, Applications of a new technique to acquire spatially resolved NMR petrophysical data, SCA Conference Paper 32, 2012. [341] W.E. Kenyon, D.F. Allen, N.V. Lisitza, Y.Q. Song, Better pore-size distributions from stimulated-echo NMR lab measurements using magnetic susceptibility contrast and small encoding angles, in: SPWLA 43rd Annual Logging Symposium, June 2–5, 2002. [342] E. Johannesen, J. Howard, A. Graue, Evaluation of wettability distributions in experimentally aged core, SCA Conference Paper Number 17, 2008. [343] E.J. Candès, J. Romberg, T. Tao, Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information, IEEE Trans. Inform. Theory 52 (2006) 489–509. [344] D. Xiao, B.J. Balcom, Two-dimensional T2 distribution mapping in rock core plugs with optimal k-space sampling, J. Magn. Reson. 220 (2012) 70–78. http://dx.doi.org/10.1016/j.jmr.2012.04.003. [345] S.G. Hyberts, A.G. Milbradt, A.B. Wagner, H. Arthanari, G. Wagner, Application of iterative soft thresholding for fast reconstruction of NMR data non-uniformly sampled with multidimensional Poisson Gap scheduling, J. Biomol. NMR 52 (2012) 315–327. [346] J.L. Paulsen, H. Cho, G. Cho, Y.Q. Song, Acceleration of multi-dimensional propagator measurements with compressed sensing, J. Magn. Reson. 213 (2011) 166–170. [347] D. Donoho, Compressed sensing, IEEE Trans. Inform. Theory 52 (2006) 1289–1306. [348] M. Lustig, D. Donoho, J.M. Pauly, Sparse MRI: the application of compressed sensing for rapid MR imaging, Magn. Reson. Med. 58 (2007) 1182–1195. [349] M. Fornasier, H. Rauhut, Iterative thresholding algorithms, Appl. Comput. Harmon. Anal. 25 (2008) 187–208. [350] S.J. Kim, K. Koh, M. Lustig, S. Boyd, D. Gorinevsky, A method for large-scale ℓ1 -regularized least squares, IEEE J. Sel. Top. Signal. 1 (2007) 606–617. [351] W. Yin, S. Osher, D. Goldfarb, J. Darbon, Bregman iterative algorithms for ℓ1 -minimization with applications to compressed sensing, SIAM J. Imaging Sci. 1 (2008) 143–168. [352] M. Lustig, Sparse MRI, Ph.D. Thesis, Stanford University, 2008. [353] G.L. Bretthorst, C. Hung, D.A. DAvignon, J.J.H. Ackerman, Bayesian analysis of time-domain magnetic resonance signals, J. Magn. Reson. 79 (1988) 369–376. [354] D. Xing, S.J. Gibbs, J.A. Derbyshire, E.J. Fordham, T.A. Carpenter, L.D. Hall, Bayesian analysis for quantitative NMR flow and diffusion imaging, J. Magn. Reson. B 106 (1995) 1–9. [355] D.J. Holland, A. Blake, A.B. Tayler, A.J. Sederman, L.F. Gladden, A Bayesian approach to characterising multi-phase flows using magnetic resonance: application to bubble flows, J. Magn. Reson. 209 (2011) 83–87. [356] T. Bayes, R. Price, An essay towards solving a problem in the doctrine of chances, Phil. Trans. 53 (1793) 370–418. [357] A. Cardenas-Blanco, C. Tejos, P. Irarrazaval, I. Cameron, Noise in magnitude magnetic resonance images, Concept. Magn. Reson. 32A (2008) 409–416. [358] G.T. Nolan, P.E. Kavanagh, The size distribution of interstices in random packings of spheres, Powder Tech. 78 (1994) 231–238. [359] N. Reboul, E. Vincens, B. Cambou, A statistical analysis of void size distribution in a simulated narrowly graded packing of spheres, Granul. Matter 10 (2008) 457–468.