Magnetic resonance for downhole complex-lithology earth formation evaluation

Many petroleum-bearing reservoirs are known as complex-lithology, complex-mineralogy earth formations where the existing nuclear magnetic resonance (NMR) analysis models require modification and extension to work properly. This paper discusses the effects of mineralogy variation in complex-lithology formation rocks that will cause NMR response variations. In particular, the existence of iron-rich authigenic clays and the nonquartz grains in siliciclastics affect the surface relaxivity and the internal field inhomogeneity. Using a simplified pore-lining clay model, we estimated that the magnitude of the internal gradient in siliciclastics is comparable with or greater than that of the instrument-generated field gradient common to the currently used NMR logging tools. To account for these effects for reservoir fluid identification and quantification from NMR data, an inversion model is created, which includes the mean susceptibility of the grain minerals, the correlation between pore size and the relaxation time of the wetting-phase fluid in the pores and the geometric restriction to the fluid molecular diffusion.


Introduction
In situ nuclear magnetic resonance (NMR) measurements in downhole environments encounter many special problems that are unseen in the laboratory and medical NMR, and pose special requirements in instrumental design as well as many limitations on the measurement techniques. In general, the in situ measurements are based on three NMR physical principles: • The proportionality of the proton NMR signal strength to the fluid volume in the voids (i.e. pore space) of the porous rock.
• The relationship of the wetting-phase fluid relaxation time to the individual pore dimensions.
• The use of relaxation times and diffusion contrasts among the fluid phases to discern reservoir fluids and non-native fluids, such as drilling fluids filtrating to the rock.
New petroleum resources are often found in complex-lithology formations. Rocks of complexlithology formations are often highly heterogeneous in texture, mineral composition and distribution, and pore and grain sizes. Physical and chemical changes can accompany the deposition of the sediments and/or post-depositional diagenetic process of the rocks. Some of these changes can result in different NMR responses to the fluids (e.g. oil, water and gas) occupying the pore spaces. Separating the NMR response due to rock geochemical and textural effects from the intrinsic fluid effects is the main emphasis today for NMR petrophysicists-to quantify saturating fluids and to assess the reservoir quality. This paper starts with a brief description of the downhole logging environments during the while-drilling and after-drilling-and-before-casing stages and their unique requirements and limitations on NMR instruments and data-acquisition methods. This part of the text is primarily for readers who are not familiar with NMR logging techniques. The discussion of the hardware and pulse sequences is focused on the tool-generated and rock mineral-grain-generated internal gradients that are the background for understanding the need for a new data inversion model. The remaining text focuses on the description and discussion of how the NMR responses in complex-lithology rocks are expected to differ from those of clean glass beads or quartzous sandstones. In this paper, we mainly limit our studies to iron-rich authigenic clays and other multimineral siliciclastic formations. We estimate the magnitude of the internal gradient in 3 B0 r B0 r NMR NMR pores in authigenic-clay-rich sandstone and take into account the duality of surface relaxivity in a feldspathic formation. To account for the internal gradient and restricted diffusion effects on quantifying the multiphase fluids in rocks, we developed a new two-dimensional (2D) NMR data inversion method that incorporates the magnetic susceptibility and the restricted diffusion. The internal gradient in the model is linked to the average susceptibility and is allowed to vary in different sized pores.

Nuclear magnetic resonance downhole logging instruments
In situ NMR earth formation measurements are mainly used during the process of drilling a well, known as logging-while-drilling (LWD) and after-drilling, but before casing the well, known as wireline (WL) openhole (OH) logging. Downhole NMR logging instruments use inside-out, antenna-type designs, so the instrument can fit in a wellbore, while the NMR-sensitive volume is in the rock formation (Coates et al 1991, Kleinberg 2001, Kruspe et al 2000. Figure 1 illustrates a centralized (left) and a side-looking (right) NMR logging tool in a borehole (Chen et al 2003), which are typical WL NMR logging tool designs. The large circular ring farthest from the NMR instrument illustrates the sensitive volume of the measurement. Colors represent radiofrequency (RF) illumination, with blue as the lowest and red as the highest B 1 strength. Both B 0 and B 1 fields are a function of the distance from the NMR logging tool, and they decrease rapidly as the distance r increases. Consequently, a typical NMR logging tool measures the formation only a few inches from the borehole. The dependence of the magnetic field on the distance r, illustrated in figure 1, limits the thickness of the slice of the NMR-sensitive volume, but allows selectivity of the depthof-investigation by choosing a proper frequency that resonates at B 0 (r ). Thus, for a long-T 1 fluid system, one can conduct multiple frequency measurements within the time required for protons to fully polarize, with each frequency irradiating protons in nonoverlapping sensitive volumes. This compensates for the poor signal-to-noise ratio due to the thin individual slice of the sensitive volume. Moreover, the static magnetic field gradient, G(r ) = − dB 0 (r ) dr , can be used to measure fluid diffusivity for discerning and quantifying multiphase fluids saturating the reservoir rock formation.
LWD NMR measurements are carried out when the drilling instruments are rotating rapidly. Therefore, unlike the WL NMR instrument, the NMR B 0 and B 1 fields for LWD instruments are desirable in azimuthally symmetric configuration. On the other hand, the penetration rate in the direction along the borehole trajectory is relatively slow compared with the logging speed in WL operations. Therefore, the prepolarization required is noncritical for LWD NMR instruments.
In addition, the drilling operation produces irregular vibrations substantially perpendicular to the borehole trajectory. If the vibration occurs in a thin-slice sensitive volume defined by a strong gradient magnetic field, each RF pulse in a sequence such as Carr-Purcell-Meiboom-Gill (CPMG) (Meiboom and Gill 1958) may irradiate substantially nonoverlapping slices of sensitive volume in a rock formation. This can cause significant, nonpredictable errors in relaxation time analysis. Therefore, a substantially uniform magnetic field is preferred in order to overcome the effect of vibrations on the NMR. Figure 2 illustrates the magnetic field configuration of the MagTrak TM instrument, an LWD NMR tool (Kruspe et al 2000). Such a design generates a relatively uniform magnetic field in the sensitive volume ring. Compared with the slice thickness of 1 mm for a typical WL NMR gradient instrument such as MReX TM , the sensitive volume slice thickness of MagTrak is 2.0 cm, with the average static magnetic field gradient being approximately 2 Gauss cm −1 in its sensitive volume.
The substantially lower magnetic field gradient of the LWD tool is insensitive to diffusion measurement for most reservoir fluids; therefore, it limits the diffusion-contrast-based fluiddiscrimination measurements. On the other hand, because the tool-generated field gradient is low, the internal gradient effect can be more observable even at a relatively low value, which can complicate the analysis of data acquired with low-gradient tools.
porosity as well as discerning multiphase fluids, measurement techniques that can maximize the relaxation time and diffusion contrasts on which the fluid discrimination and irreducible porosity are based have been the focus of NMR logging pulse sequence designs. Among these, a variation of CPMG sequence with assorted wait time TW, interecho time TE and the echo train length NE is most widely used with the high-gradient, WL OH NMR logging instruments. This allows multiple-frequency operation and the applied static field gradient decreases with frequency, which makes G(f ) the additional controllable variable. A typical multiple echo train sequence for the NMR WL OH logging tool can be found in the literature on NMR logging (e.g. Freedman andHeaton 2004, Chen et al 2009). The data acquired by such a sequence are expressed as where M is the echo amplitude, a is the partial porosity on the grid with index (i, j, k), D is the fluid diffusivity, γ is the gyromagnetic ratio, G is the combination of the magnetic field gradient generated by the NMR logging instrument and the internal field gradient and R = T 1 /T 2 . The indices I, J and K represent the numbers of grids along the T 2 , D and R dimensions, respectively. Multiple echo trains are often processed together to obtain multidimensional parameter-space results (D, T 2 , R) from which the fluid quantification is determined , 2009, DiRosa et al 2006, Heaton et al 2008, Hursan and Chen 2005, Song et al 2002, Sun et al 2006. For a low-gradient LWD NMR instrument such as MagTrak, fluid discrimination is largely based on T 1 or T 2 contrasts. To make diffusion contrast sensitive, one must use very large TE. This is impractical even for gas, as the fast-relaxing components, usually representing porosity associated with clay-bound water, will be lost by a large TE. To overcome such a problem, Chen (2009) proposed an echo-decay-acceleration method utilizing a pulse sequence illustrated in figure 3, where two echo trains are included. The first echo train is the normal CPMG echo train with a short TE = TE S . The second echo train is acquired with two TE values. The first n S echoes are acquired with the short TE S the same as Echo Train 1, where a typical length of the 6 first portion of the second echo train, t S = n S · TE S , is about 100 ms. We will go back to the t S later but for now we will complete the echo train description first. After t S the TE value is switched to a long TE value (TE L ). The TE L value is selected such that the combination of G · TE L is sufficiently sensitive to distinguish gas from liquids based on their diffusion contrast. The second portion of the echo train has the length t L = n L · TE L and typically ranges from several hundreds of ms to 1 s. As a point of clarification, the two portions of the second echo train form a single, continuous decay curve.
The initial part of the echo train decay, acquired with TE S , is represented with the highest sampling rate, so that the fast decay components, normally associated with bound fluids (often called BVI in the oil industry), can be reliably determined. The slowly decaying gas component, on the other hand, experiences little reduction in amplitude at t = t S .
For example, for 100% methane at 75 • C and 3000 psi, the gas T 2intr ≈ T 1 = 2734 ms. Using a typical MagTrak field gradient value of 2 Gauss cm −1 and TE S = 0.6 ms, the decay of gas signal at t S = 100 ms is less than 5%, calculated using and In contrast, for a typical irreducible water relaxation time of 30 ms, only 3.5% BVI signal is left at 100 ms. Even if the t S is doubled to 200 ms, 91% of the gas signal should remain at the end of the first portion of the echo train, i.e. t = 200 ms. The above calculation has assumed zero internal gradient. When the internal gradient is comparable with or stronger than the toolgenerated gradient, the result can be significantly different, potentially causing misidentification of reservoir fluids.

Mineralogy effect on relaxation time distribution in complex lithology formations
NMR relaxation time (T 1 or T 2 ) distribution data are widely used as the entry point for many NMR-based formation property interpretation models. The rationale is that the relaxation rate T −1 1 of a wetting-phase fluid in a single pore is considered to be proportional to the product of the surface-to-volume ratio (S/V ) and the surface relaxivity (ρ), where T 1bulk is the bulk relaxation time of the fluid and the index i represents the ith pore. A similar expression exists for T 2 . Strictly speaking, the surface relaxivity for T 1 and T 2 may be slightly different. The pore-surface-to-pore-volume ratio is equivalent to the inverse of the average curvature of the pore, which is defined as the characteristic pore size. Reservoir rocks often have a distribution of pore sizes; thus a distribution of T 1 and T 2 is expected. When the wetting phase in the formation rock is water, T 1 and T 2 of water are in the order of several seconds at reservoir temperatures; thus the second term in (4) dominates.
The surface relaxivity, ρ, plays the role of the proportionality constant between wettingphase fluid T 2 and the pore size. However, as ρ depends on the fluid and rock matrix minerals, it is not a constant for complete lithology formation containing multimineral grains and/or for rock and fluids that are not strongly water-wet. For a multimineral formation where each pore is surrounded by grains of the same mineral but are all of the same wetting characteristics, (4) should be more correctly written as if the ith pore is bordered by grains of the kth-type mineral. On the other hand, if the ith pore is bordered by grains of n different types of minerals, the above equation is modified to where Clearly, one needs mineralogy information to correlate the nonlinear behavior between T 2 distribution and pore-size distribution in complex-lithology situations. Figure 4 shows two thin-section photographs of two rock samples taken from two different sections of a Baker Hughes beta test well located in Oklahoma. The left photograph shows a quartzous sandstone formation and the right feldspathic sandstone formation. The corresponding T 2 distributions of these two sections are shown in figure 5. The position of the long T 2 relaxation time peak shift in this figure can be interpreted from the pore-size difference for pores surrounded by quartzous grains. The fast-relaxing T 2 component in the feldspathic section is contributed by fluids in smaller pores surrounded by feldspathic grains. As ρ feldspar and ρ quartz might be different, an adjustment of the relaxivity value is needed in order to accurately translate the bimodal T 2 distribution to pore-size distribution. In this case, ρ feldspar /ρ quartz = 3.3 is used, with the assumption that the same surface roughness factor applies to both mineral grains. Had we used a single ρ quarz to compute pore-size distribution, the pore sizes associated with the pores surrounded by the feldspar grains would have been overestimated by a factor of 3.3. Knowledge of the mineral composition enables consideration of these nonlinear effects.
For authigenic clay-lining pores, a single surface relaxivity can be used in the inversion model. For siliciclastic rocks containing multimineral grains such as the quartzous-feldspathic rocks, quantitative volumes of the mineral grains are needed for assigning the corresponding surface relaxivity values in the fractions of pore spaces.

Effect of authigenic clay on the internal gradient
Many complex-lithology formations are clay-rich clastic formations. The internal field gradient in pores with iron-rich authigenic clays can be greater than the logging tool-generated field gradient strength. Authigenic clays are formed in pores from minerals existing in the connate water. The growth process often results in clay spikes or flakes forming on the surface of the pores; the pore-lining chlorite flake has a much higher magnetic susceptibility than the surrounding pore fluids. The geometry and distribution of these clays are highly variable. It is practically impossible and unnecessary to track the individual clay flakes to calculate the magnetic field inhomogeneity. Instead, one can use a simplified geometry to simulate the internal gradient and, subsequently, estimate the average effect on NMR measurements. Zhang et al (2001) approximates the pore-lining clay flakes to rectangular clay flakes, with the clay/fluid interfaces either parallel or perpendicular to the homogeneous magnetic field B 0 . They used Green's function representation to calculate the internal gradient magnitude for individual clay flakes having a magnetic susceptibility χ. Then, the field gradient caused by multiple clay flakes is computed by superposition of Green's function. The use of a square pore system simplifies the 2D model and yet is sufficient for modeling clays. A ratio of 7 : 1 between the length and width of a typical chlorite clay flake is chosen and the gap between two clay flakes is set as the clay width. The proportion of the macropores relative to the micropores is modeled by the number of clay flakes on each side of the square pore. Using a hypothetical 0.32 irreducible water saturation for computing fraction of micropore volume, 15 pore-lining clay flakes are included in the model. Figure 6 shows the contour lines of dimensionless gradients for the micropore and macropore of this model authigenic clay lined pore. The values of the contour lines differ by a factor of two. High gradients occur around the clay tips independent of the orientations of Also shown in figure 6 are the minimum, maximum, mean and standard deviation values of the dimensionless gradient. Dimensional gradient values in units of Gauss cm −1 can be easily determined from dimensionless values through the width of clay flakes. Using a clay width of 0.2 µm, the contour lines of 0.01 and 0.16 are approximately 2 and 32 Gauss cm −1 , respectively, compared with the field gradients of 15-40 Gauss cm −1 of a typical WL logging tool and 2 Gauss cm −1 of the MagTrak LWD tool. Clearly, the magnetic field within North Burbank (NB) sandstone is not homogeneous; the use of a single 'average' gradient value to describe such a complex pore system will be a simplistic approach. Zhang (2001) has compared the TE-dependent T 2 relaxation times in NB sandstones with the average internal gradient estimates from the simplified model. It has been concluded that the estimated internal gradient using the simplified model is consistent with the experimental measurements to the correct order-of-magnitude. However, no linear (TE/2) 2 dependence is observed. While it is neither necessary nor possible to account for the actual pore and clay geometries in the model, one must also take into account the fact that the average gradient is dependent on pore dimension, which for the wetting-phase saturated pores is also proportional to the NMR relaxation time of the fluid in the pore. Moreover, the geometric restriction to the molecule diffusion needs to be taken into account. In the next section, we describe a new approach that includes all these effects (the magnetic susceptibility variation, pore-sizedependent internal gradient and geometric restriction to the diffusion) in the NMR data inversion model.

Accounting for internal gradient and restricted diffusion for two-dimensional (2D) NMR inversion
Internal gradient effects coupled with restricted diffusion have been well-recognized but longstanding problems in NMR data inversion. To quantify the influence of the internal gradient and restricted diffusion effects in formation and multiphase fluids discernment, we start with the Bloch-Torrey equation (Torrey 1956) which describes the magnetization evolution, m, under a field gradient, G. Under pore geometry confinement, a proper boundary condition, will also account for the restricted diffusion effect. The complexity of such effects arises from the fact that the equation is not exactly solvable, unless under extreme cases such as constant gradient with no boundaries (free diffusion). Numerically, however, and with simple boundary geometries and a known gradient profile, the Block-Torrey equation can be solved by the finitedifference method (Blees 1994) or the multiple correlation function (MCF) approach developed by Grebenkov (2007). Alternatively, the Monte Carlo technique can be used to simulate the Brownian motion of individual spins and calculate the accumulated evolution of magnetization. Our approach takes advantage of the discrete nature of the inversion algorithm and calculates the responses under every possible condition (every combination of T 2 and D, with assigned magnetic susceptibility and surface relaxivity). The relaxation time rate in this case is a function of various variables, where B m , G m are the magnetic field and field gradient generated by the logging instrument, TE n represents one instance of TE as configured in the pulse sequence, ρ k and χ l are possible numbers for surface relaxivity and magnetic susceptibility, respectively, and D i and T 2 j represent one point in the 2D inversion grid. As all variables take discrete numbers, the possible values of R become a very large matrix, although it is still manageable and can be numerically precomputed and stored in memory. At the time of inversion, the values are fetched from memory efficiently. For details of this method, see the work by Li and Chen (2010).
Owing to the heterogeneous nature of the rock and the lack of a high-resolution and representative image of the rock, it is practically impossible to model the actual pore geometry of the rock. Consequently, some assumptions and simplifications are made. A 1D slab geometry and a 3D spherical geometry were both used for numerical calculation and the 3D model was eventually chosen as a statistical model for the internal gradient profile in a rock formation. When diffusion is restricted between parallel planes (1D slab), we have used the following model for the internal field: where L is the length between the two parallel planes. In the case of spherical boundaries (3D sphere), the following internal field distribution along the z-direction was used: The above model satisfies the following conditions: (i) the maximum and minimum of the internal field are proportional to χ B 0 , where χ is the magnetic susceptibility contrast; (ii) the field gradient is generally higher at the liquid-solid interface; (iii) the profile of the internal field does not change with the size of the pores, or in other words, the magnitude of the internal field gradient is inversely proportional to the pore size. We assume that the mineralogy distribution in the rock can be treated as uniform for the purpose of this internal field gradient estimation. We use the MCF approach to calculate the behavior of our model system. For details of the calculation, see the review paper by Grebenkov (2007). To verify the validity of our approach, we also conducted a random-walk simulation in a pore system of randomly packed spheres (as grains), with the magnetic field generated by 10 000 magnetic dipoles (grains). The magnetic dipole moment of each sphere is where R is the radius of the sphere, χ is the difference in magnetic susceptibility between the grain material and the fluid, B 0 is the external magnetic field and µ 0 is the permeability of free space. The magnetic field at position r due to the contribution of the ith sphere is Figure 8. Inversion results using logging data. Left: using the original algorithm; right: using the algorithm that incorporates the internal gradient and restricted diffusion. Improvements are dramatic, and we believe that the graph on the right also gives a more accurate quantification of the gas signal, as the water signal appears on the water diffusivity line and the gas and water signals are separated clearly.
where δr i = r − r i is the vector from the sphere position r i to r. If we only consider B(r) in the z-direction and assume that B 0 lies along the z-direction, Summing B i z (r) over 10 000 spheres surrounding the pore will help determine the internal field at any position r in the pore.
Details of Monte Carlo simulation can be found in the work by Valckenborg et al (2002). Figure 7 shows a comparison of the results of Monte Carlo simulation and those of numerical calculation based on a simplified model. We can see that our model system statistically captures the essential features of the internal gradient and demonstrates similar relaxation behavior as that of a more realistic pore system, with more realistic internal field distribution. The advantage of our model is that it enables tremendous time savings for the calculation of each T 2 , D, χ and ρ (surface relaxivity) combination. Of course the distribution of the internal gradient in a real rock pore system is more complicated than that in a random spherical pack of magnetic dipoles. We believe that as long as the length scale of the internal gradient variation is the same as the length scale of the pore, our simplified spherical model can be applied. Examples of such cases can be the case when the paramagnetic materials on the grain surface are smoothly distributed or are sufficiently close to one another or the case when each pore only has one or two paramagnetic centers. When the pore has several paramagnetic centers and the length scale of the internal gradient is significantly smaller than the pore size, the inversion results from our simplified model may not be that satisfactory any longer. Further study is planned on this subject. Figure 8 shows a comparison of different processing methods using test well data. As a general rule, in the inversion process we choose a magnetic susceptibility contrast such that the water signal lies roughly along the water diffusivity line. The surface relaxivity of 10 µm s −1 is 13 normally used for sandstone and this can be changed. We have found that the surface relaxivity does not affect the result very much unless a drastic change of more than one order of magnitude is made. A different susceptibility contrast may change the relative position of the spots but has less impact on the quantification result.
The improvement shown in figure 8 is significant. The data are acquired from a complexlithology gas-bearing formation interval. Without inclusion of the internal gradient and restricted diffusion effects in the inversion process (left) the water signal is not on the waterdiffusivity line (showing as a thin white horizontal line). After taking into account these effects, the inversion yields a result that is much more reasonable for a water-bearing reservoir.

Summary and conclusions
We have discussed the different logging environments in LWD and WL stages of reservoir exploration and the impact on NMR measurements and hardware preferences and constraints in these environments. Particular attention was given to the comparison between the tool-generated field gradients and the internal field gradients. The effects of multiple minerals in complexlithology formations are described in terms of surface relaxivity and internal gradients. An inversion model incorporating internal gradient and restricted diffusion is developed. In this model, the internal gradient is allowed to vary due to pore size and is linked to the relaxation times of the wetting-phase pore fluid. The effect of magnetic susceptibility variation due to multimineral grains is treated by an average value. Using this approach, we are able to improve reservoir multiphase fluid discrimination, which is particularly useful in complex-lithology reservoirs.