Methods and systems for determining surface relaxivity of a medium using nuclear magnetic resonance

Abstract

Methods and systems for determining surface relaxivity from nuclear magnetic resonance measurements relate to applying multiple nuclear magnetic resonance (NMR) diffusion editing Carr-Purcell-Meiboom-Gill (CPMG) pulse sequences to the porous medium, wherein the diffusion editing CPMG pulse sequences have a diffusion encoding time Δ; receiving NMR data generated by the pulse sequences; processing the received NMR data to obtain a distribution f(T.sub.2,D) for the diffusion encoding time Δ; repeating the applying, the receiving, and the processing at least one time for pulse sequences having different respective diffusion encoding times Δ to obtain respective distributions f(T.sub.2,D) corresponding respectively to the different diffusion encoding times Δ; and utilizing the respectively obtained distributions f(T.sub.2,D) to generate a surface relaxivity (ρ) determination.

Claims

1. A method of determining surface relaxivity of a porous medium, comprising: applying multiple nuclear magnetic resonance (NMR) diffusion editing Carr-Purcell-Meiboom-Gill (CPMG) pulse sequences to the porous medium, wherein the diffusion editing CPMG pulse sequences have gradient pulses for diffusion editing separated by a diffusion encoding time Δ; receiving NMR data generated by the pulse sequences; processing the received NMR data to obtain a distribution f(T.sub.2,D) for the diffusion encoding time Δ, where T.sub.2 is transverse relaxation time and D is diffusion coefficient; repeating the applying, the receiving, and the processing at least one time for pulse sequences having different respective diffusion encoding times Δ to obtain respective distributions f(T.sub.2,D) corresponding respectively to the different diffusion encoding times Δ; and utilizing the respectively obtained distributions f(T.sub.2,D) together to generate a surface relaxivity (ρ) determination.

2. The method according to claim 1, further comprising: calculating a mean diffusion coefficient D(T.sub.2s) for each T.sub.2 value of each distribution f(T.sub.2,D), where T.sub.2s is the transverse relaxation time due to surface relaxation; and calculating a respective T.sub.2 distribution g(T.sub.2s) for each distribution f(T.sub.2,D), wherein the surface relaxivity (ρ) determination is generated utilizing the calculated T.sub.2 distributions and the calculated mean diffusion coefficients (T.sub.2s).

3. The method according to claim 1, further comprising: calculating a weighted mean transverse relaxation time T.sub.2s(D) for each D value of each distribution f(T.sub.2,D) to compute a functional relation between D and T.sub.2, wherein the surface relaxivity (ρ) determination is generated utilizing the functional relation between D and T.sub.2.

4. The method according to claim 1, wherein the generating multiple nuclear magnetic resonance CPMG pulse sequences comprises determining the diffusion encoding time Δ such that the characteristic diffusion length of a saturated fluid expected to be disposed in the porous medium corresponds to an expected pore size of the porous medium.

5. The method according to claim 1, wherein: the utilizing comprises conducting a data-fitting error minimization.

6. The method according to claim 5, wherein: the data-fitting error minimization comprises choosing values for the surface relaxivity and at least one parameter to cause a fitting function of the surface relaxivity and the at least one parameter to closely approximate a data set of the average diffusion coefficients D(Δ,T.sub.2s).

7. The method according to claim 6, wherein: the data-fitting error minimization comprises utilizing a logarithmic error.

8. The method according to claim 7, wherein: the logarithmic error is defined according to err=Σ.sub.ΔT.sub.2s{g(Δ,T.sub.2s)*[log D.sub.fit(Δ,T.sub.2s)−log D(Δ,T.sub.2s)]}.sup.2, where D.sub.fit(Δ,T.sub.2s) is the fitting function.

9. The method according to claim 8, wherein: the fitting function is defined by $\frac{D\left(\Delta ,T_{2s}\right)}{D_0}=1-\left(1-\alpha \right)\frac{\beta \frac{L_D}{\rho T_{2s}}+\left(1-\alpha \right){\left(\frac{L_D}{\mathrm{\Gamma \rho }{T}_{2s}}\right)}^2}{\left(1-\alpha \right)+\beta \frac{L_D}{\rho T_{2s}}+\left(1-\alpha \right){\left(\frac{L_D}{\mathrm{\Gamma \rho }{T}_{2s}}\right)}^2}$ where D.sub.0 is a bulk diffusion coefficient, $\alpha =D\left(\infty \right)/D_0,\beta =\frac{4}{9\sqrt{\pi }},L_D=\sqrt{D_0\Delta },\mathrm{and}\Gamma =L_M/\rho T_{2s},$ where L.sub.M is a length scale characterizing a transition between a short-time and long-time diffusion limit.

10. The method according to claim 9, wherein: the at least one parameter comprises D.sub.0, α, and Γ.

11. The method according to claim 1, wherein: the processing the resulting data comprises solving $M\left(q_i,{\tau }_j\right)={\Sigma }_{k,l}f\left(T_{2,k},D_l\right)e^{-q_i^2{D}_l\Delta }{e}^{-\frac{{\tau }_i}{T_{2,k}}}$ to obtain the distribution f(T.sub.2,D) where M(q.sub.i,τ.sub.j) is the resulting data, q is the encoding strength, and τ is the time in CPMG decay.

12. The method according to claim 11, wherein: the solving comprises using a fast Laplace inversion.

13. The method according to claim 11, wherein: the calculating a mean diffusion coefficient D(T.sub.2s) for each T.sub.2 value comprises converting T.sub.2 indices into T.sub.2s indices by relabeling the T.sub.2 index as a T.sub.2s index utilizing $\frac{1}{T_{2s}}=\frac{1}{T_2}-\frac{1}{T_{2b}}$ where T.sub.2b is the T.sub.2 value of bulk fluid in the porous medium.

14. The method according to claim 13, wherein: the calculating a mean diffusion coefficient D(T.sub.2s) for each T.sub.2 value further comprises calculating according to $\stackrel{\_}{D}\left(\Delta ,T_{2s}\right)=\frac{{\Sigma }_{i}f\left(T_{2s},D_i\right)*D_i^2}{{\Sigma }_{i}f\left(T_{2s},D_i\right)*D_i}.$

15. The method according to claim 11, wherein: the calculating a mean diffusion coefficient D(T.sub.2s) for each T.sub.2 value further comprises calculating according to $\stackrel{\_}{D}\left(\Delta ,T_{2s}\right)=\sqrt{\frac{{\Sigma }_{i}f\left(T_{2s},D_i\right)*D_i^2}{{\Sigma }_{i}f\left(T_{2s},D_i\right)}}.$

16. The method according to claim 11, wherein: the calculating a mean diffusion coefficient D(T.sub.2s) for each T.sub.2 value further comprises calculating according to $\stackrel{\_}{D}\left(\Delta ,T_{2s}\right)=\frac{{\Sigma }_{i}f\left(T_{2s},D_i\right)*D_i}{{\Sigma }_{i}f\left(T_{2s},D_i\right)}.$

17. The method according to claim 11, wherein: the calculating a T.sub.2 distribution g(T.sub.2s) comprises calculating according to g(Δ,T.sub.2s)=Σ.sub.D.sub.if.sub.Δ(T.sub.2s,D.sub.i).

18. The method according to claim 9, wherein: the processing the resulting data comprises solving $M\left(q_i,{\tau }_j\right)={\Sigma }_{k,l}f\left(T_{2,k},D_l\right)e^{-q_i^2{D}_l\Delta }{e}^{-\frac{{\tau }_i}{T_{2,k}}}$ using a fast Laplace inversion to obtain the distribution f(T.sub.2,D) where M(q.sub.i,τ.sub.j) is the resulting data, q is the encoding strength, and τ is the time in CPMG decay.

19. The method according to claim 9, wherein: the calculating a mean diffusion coefficient D(T.sub.2s) for each T.sub.2 value comprises converting T.sub.2 indices into T.sub.2s indices by relabeling the T.sub.2 index as a T.sub.2s index utilizing $\frac{1}{T_{2s}}=\frac{1}{T_2}-\frac{1}{T_{2b}}$ where T.sub.2b is the T.sub.2 value of bulk fluid in the porous medium.

20. The method according to claim 19, wherein: the calculating a mean diffusion coefficient D(T.sub.2s) for each T.sub.2 value further comprises calculating according to $\stackrel{\_}{D}\left(\Delta ,T_{2s}\right)=\frac{{\Sigma }_{i}f\left(T_{2s},D_i\right)*D_i^2}{{\Sigma }_{i}f\left(T_{2s},D_i\right)*D_i}.$

21. The method according to claim 19, wherein: the calculating a mean diffusion coefficient D(T.sub.2s) for each T.sub.2 value further comprises calculating according to $\stackrel{\_}{D}\left(\Delta ,T_{2s}\right)=\sqrt{\frac{{\Sigma }_{i}f\left(T_{2s},D_i\right)*D_i^2}{{\Sigma }_{i}f\left(T_{2s},D_i\right)}}.$

22. The method according to claim 19, wherein: the calculating a mean diffusion coefficient D(T.sub.2s) for each T.sub.2 value further comprises calculating according to $\stackrel{\_}{D}\left(\Delta ,T_{2s}\right)=\frac{{\Sigma }_{i}f\left(T_{2s},D_i\right)*D_i}{{\Sigma }_{i}f\left(T_{2s},D_i\right)}.$

23. The method according to claim 19, wherein: the calculating a T.sub.2 distribution g(T.sub.2s) comprises calculating according to g(Δ,T.sub.2s)=Σ.sub.D.sub.if.sub.Δ(T.sub.2s,D.sub.i).

24. The method according to claim 1, wherein: the repeating a)-d) multiple respective times comprises repeating a) multiple times and storing the acquired resulting data prior to the processing, the calculating a mean diffusion coefficient, and the calculating a T.sub.2 distribution.

25. The method according to claim 1, further comprising: displaying an indication of the determination of surface relaxivity (p).

26. The method according to claim 25, wherein: the displaying comprises generating a log of surface relaxivity as a function of depth in a borehole.

27. A method of determining surface relaxivity of a porous medium, comprising: (a) generating multiple nuclear magnetic resonance (NMR) diffusion editing—CPMG (Carr-Purcell-Meiboom-Gill) pulse sequences with a diffusion encoding time Δ that interact with the porous medium, and acquiring resulting data; (b) processing, by a processor, the resulting data to obtain a distribution f(T.sub.2,D) for that Δ; (c) calculating, by a processor, a mean diffusion coefficient D(T.sub.2s) for each T.sub.2 value, where T.sub.2s is the transverse relaxation time due to surface relaxation; (d) calculating, by a processor, a T.sub.2 distribution g(T.sub.2s); and (e) utilizing, by a processor, the calculated T.sub.2 distribution and the calculated mean diffusion coefficient in order to generate a determination of surface relaxivity (ρ) by choosing values for the surface relaxivity and at least one parameter to cause a fitting function of the surface relaxivity and a parameter Γ to approximate the average diffusion coefficient D(T.sub.2s), the fitting function is defined by $\frac{D\left(\Delta ,T_{2s}\right)}{D_0}=1-\left(1-\alpha \right)\frac{\beta \frac{L_D}{\rho T_{2s}}+\left(1-\alpha \right){\left(\frac{L_D}{\mathrm{\Gamma \rho }{T}_{2s}}\right)}^2}{\left(1-\alpha \right)+\beta \frac{L_D}{\rho T_{2s}}+\left(1-\alpha \right){\left(\frac{L_D}{\mathrm{\Gamma \rho }{T}_{2s}}\right)}^2}$ where D.sub.0 is a bulk diffusion coefficient, $\alpha =D\left(\infty \right)/D_0,\beta =\frac{4}{9\sqrt{\pi }},L_D=\sqrt{D_0\Delta },\mathrm{and}\Gamma =L_M/\rho T_{2s},$ where L.sub.M is a length scale characterizing a transition between a short-time and long-time diffusion limit.

28. The method according to claim 27, wherein: the processing the resulting data comprises solving $M\left(q_i,{\tau }_j\right)={\Sigma }_{k,l}f\left(T_{2,k},D_l\right)e^{-q_i^2{D}_l\Delta }{e}^{-\frac{{\tau }_i}{T_{2,k}}}$ to obtain the distribution f(T.sub.2,D) where M(q.sub.i,τ.sub.j) is the resulting data, q is the encoding strength, and τ is the time in CPMG decay.

29. The method according to claim 28, wherein: the calculating a T.sub.2 distribution g(T.sub.2s) comprises calculating according to g(Δ,T.sub.2s)=Σ.sub.D.sub.if.sub.Δ(T.sub.2s,D.sub.i).

30. The method according to claim 27, further comprising: displaying an indication of the determination of surface relaxivity (ρ).

31. A method utilizing nuclear magnetic resonance (NMR) measurements to determine surface relaxivity of a porous medium, comprising: (a) generating multiple NMR diffusion editing pulse sequences for each of a plurality of diffusion encoding times Δ that interact with the porous medium, and acquiring resulting data; (b) processing, by a processor, the resulting data to obtain mean diffusion coefficient D(T.sub.2s) for each T.sub.2 value of each of the plurality of encoding times Δ, where T.sub.2s is the transverse relaxation time due to surface relaxation, and to obtain T.sub.2 distributions g(T.sub.2s); and (c) utilizing, by a processor, respective calculated T.sub.2 distributions and respective calculated mean diffusion coefficients in a data-fitting error minimization procedure in order to generate a determination of surface relaxivity (ρ), where the data-fitting error minimization procedure utilizes a fitting function with a parameter that permits the length scale characterizing the transition between the short-time and long-time diffusion limit (L.sub.M) to scale with pore size.

32. The method according to claim 31, wherein: the data-fitting error minimization comprises choosing values for the surface relaxivity and at least one parameter to cause a fitting function of the surface relaxivity and the at least one parameter to closely approximate a data set of the average diffusion coefficients D(Δ,T.sub.2s).

33. The method according to claim 32, wherein: the data-fitting error minimization comprises utilizing a logarithmic error.

34. The method according to claim 33, wherein: the logarithmic error is defined according to err=Σ.sub.Δ,T.sub.2s{g(Δ,T.sub.2s)*[log D.sub.fit(Δ,T.sub.2s)−log D(Δ,T.sub.2s)]}.sup.2, where D.sub.fit(Δ,T.sub.2s) is the fitting function.

35. The method according to claim 31, wherein: the fitting function is defined by $\frac{D\left(\Delta ,T_{2s}\right)}{D_0}=1-\left(1-\alpha \right)\frac{\beta \frac{L_D}{\rho T_{2s}}+\left(1-\alpha \right){\left(\frac{L_D}{\mathrm{\Gamma \rho }{T}_{2s}}\right)}^2}{\left(1-\alpha \right)+\beta \frac{L_D}{\rho T_{2s}}+\left(1-\alpha \right){\left(\frac{L_D}{\mathrm{\Gamma \rho }{T}_{2s}}\right)}^2}$ where D.sub.0 is a bulk diffusion coefficient, $\alpha =D\left(\infty \right)/D_0,\beta =\frac{4}{9\sqrt{\pi }},L_D=\sqrt{D_0\Delta },\mathrm{and}\Gamma =L_M/\rho T_{2s},$ where L.sub.M is a length scale characterizing a transition between a short-time and long-time diffusion limit.

36. The method according to claim 35, wherein: the at least one parameter comprises D.sub.0, α, and Γ.

37. The method according to claim 31, wherein: the processing the resulting data comprises solving $M\left(q_i,{\tau }_j\right)={\Sigma }_{k,l}f\left(T_{2,k},D_l\right)e^{-q_i^2{D}_l\Delta }{e}^{-\frac{{\tau }_i}{T_{2,k}}}$ to obtain the distribution f(T.sub.2,D) where M(q.sub.i,τ.sub.j) is the resulting data, q is the encoding strength, and τ is the time in CPMG decay.

38. The method according to claim 37, wherein: the solving comprises using a fast Laplace inversion.

39. The method according to claim 31, further comprising: displaying an indication of the determination of surface relaxivity (ρ).

40. The method according to claim 39, wherein: the displaying comprises generating a log of surface relaxivity as a function of depth in a borehole.

41. A method of determining surface relaxivity of a porous medium, comprising: generating multiple nuclear magnetic resonance (NMR) diffusion editing pulse sequences having gradient pulses for diffusion editing separated for each of a plurality of diffusion encoding times Δ, wherein said multiple NMR diffusion editing pulse sequences that interact with the porous medium, and acquiring resulting data; and without utilizing non-NMR measurements, using a processor to fit a function relating the surface relaxivity to NMR diffusion coefficients and diffusion lengths to diffusion coefficients calculated from the resulting data in order to obtain a surface relaxivity determination for the porous medium.

42. The method according to claim 41, wherein: the function is defined by $\frac{D\left(\Delta ,T_{2s}\right)}{D_0}=1-\left(1-\alpha \right)\frac{\beta \frac{L_D}{\rho T_{2s}}+\left(1-\alpha \right){\left(\frac{L_D}{\mathrm{\Gamma \rho }{T}_{2s}}\right)}^2}{\left(1-\alpha \right)+\beta \frac{L_D}{\rho T_{2s}}+\left(1-\alpha \right){\left(\frac{L_D}{\mathrm{\Gamma \rho }{T}_{2s}}\right)}^2}$ where ρ is the surface relaxivity, D.sub.0 is a bulk diffusion coefficient, $\alpha =D\left(\infty \right)/D_0,\beta =\frac{4}{9\sqrt{\pi }},L_D=\sqrt{D_0\Delta }$ is a diffusion length, and Γ=L.sub.M/ρT.sub.2s, where L.sub.M is a length scale characterizing a transition between a short-time and long-time diffusion limit.

43. The method according to claim 41, further comprising: displaying an indication of the determination of surface relaxivity (ρ).

44. A nuclear magnetic resonance (NMR) system comprising: a coil for applying a NMR pulse sequence to a substance; a NMR transmitter coupled to the coil; a processor; and a memory storing instructions executable by the processor to apply multiple nuclear magnetic resonance (NMR) diffusion editing Carr-Purcell-Meiboom-Gill (CPMG) pulse sequences to the porous medium, wherein the diffusion editing CPMG pulse sequences have gradient pulses for diffusion editing separated by a diffusion encoding time Δ; receive NMR data generated by the pulse sequences; process the received NMR data to obtain a distribution f(T.sub.2,D) for the diffusion encoding time Δ; repeat the applying, the receiving, and the processing at least one time for pulse sequences having different respective diffusion encoding times Δ to obtain respective distributions f(T.sub.2,D) corresponding respectively to the different diffusion encoding times Δ; and utilize the respectively obtained distributions f(T.sub.2,D) together to generate a surface relaxivity (ρ) determination.

45. The NMR system of claim 44, wherein the NMR system is part of a wellbore tool for investigating formations.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1 is a flow chart of a method for determining the surface relaxivity of a medium.

(2) FIG. 2 is a diagram of an NMR pulse diagram sequence used for determining the surface relaxivity of a medium.

(3) FIGS. 3a-3e are DT.sub.2 distribution plots and fitting results for five different diffusion encoding times for a medium with elements having a first diameter size.

(4) FIG. 3f is a T.sub.2 distribution for the medium having a first diameter size.

(5) FIGS. 4a-4e are DT.sub.2 distribution plots and fitting results for five different diffusion encoding times for a medium with elements having a second diameter size.

(6) FIG. 4f is a T.sub.2 distribution for the medium having a second diameter size.

(7) FIGS. 5a-5e are DT.sub.2 distribution plots and fitting results for five different diffusion encoding times for a Berea sandstone rock core sample.

(8) FIG. 5f is a T.sub.2 distribution for the Berea sandstone rock core sample.

(9) FIGS. 6a-6e are DT.sub.2 distribution plots and fitting results for five different diffusion encoding times for a Fontaineblue sandstone rock core sample.

(10) FIG. 6f is a T.sub.2 distribution for the Fountainblue sandstone rock core sample.

(11) FIG. 7A shows a wireline NMR instrument deployed in a wellbore.

(12) FIG. 7B shows a logging while drilling NMR instrument deployed in a wellbore.

DETAILED DESCRIPTION

(13) Before turning to the drawings, an understanding is useful of some of the physics underlying hereinafter-described embodiments.

(14) The time dependent diffusion coefficient in porous media depends on surface volume ratio S/V. See, e.g., Mitra, P. P. et al., “Diffusion propagator as a probe of the structure of porous media,” Physical Review Letters, Vol. 68, pp. 3555-3558 (1992); and Mitra, P. P. et al., “Short-time behavior of the diffusion coefficient as a geometrical probe of porous media,” Physical Review B, Vol. 47, p. 8565 (1993). In the short time limit, i.e. where the characteristic diffusion length (√{square root over (2D.sub.0Δ)}) is much smaller than the pore size, the apparent diffusion coefficient decrease is proportional to S/V, described by

(15) $\begin{array}{cc}\frac{D\left(\Delta \right)}{D_0}=1-\frac{4}{9\sqrt{\pi }}\frac{S}{V}\sqrt{D_0\Delta },& \left(1\right)\end{array}$
where D.sub.0 is the bulk diffusion coefficient and Δ is the diffusion time. In the long time limit, D(Δ) approaches the asymptotic value D(∞)/D.sub.0=α. A Padé approximation can be used to connect the two limits,

(16) $\begin{array}{cc}\frac{D\left(\Delta \right)}{D_0}=1-\left(1-\alpha \right)\frac{\beta L_D\left(\frac{S}{V}\right)+\left(1-\alpha \right){\left(\frac{L_D}{L_M}\right)}^2}{\left(1-\alpha \right)+\beta L_D\left(\frac{S}{V}\right)+\left(1-\alpha \right){\left(\frac{L_D}{L_M}\right)}^2}& \left(2\right)\end{array}$
Where

(17) $\beta =\frac{4}{9\sqrt{\pi }},L_D=\sqrt{D_0\Delta },$
and L.sub.M is the length scale characterizing the transition between the short-time and long-time diffusion limit. See, Latour, L. L., et al., “Time-dependent diffusion coefficient of fluids in porous media as a probe of surface-to-volume ratio,” Journal of Magnetic Resonance, Series A, Vol. 101 pp. 342-346 (1993); and Hürlimann, M., et al., “Diffusion-relaxation distribution functions of sedimentary rocks in different saturation states,” Magnetic Resonance Imaging, Vol. 21, pp. 305-310 (2003).

(18) NMR T.sub.2 relaxation time is also affected in porous media by the surface to volume ratio as described by the expression

(19) $\begin{array}{cc}\frac{1}{T_{2s}}=\rho \frac{S}{V}& \left(3\right)\end{array}$
where T.sub.2s is the transverse relaxation time due to surface relaxation (where

(20) $\frac{1}{T_{2s}}=\frac{1}{T_2}-\frac{1}{T_{2b}},$
with T.sub.2b being the bulk relaxation) and ρ is the surface relaxivity. Combining equations (2) and (3) leads to a formula relating T.sub.2s and D(Δ):

(21) $\begin{array}{cc}\frac{D\left(\Delta \right)}{D_0}=1-\left(1-\alpha \right)\frac{\beta \frac{L_D}{\rho T_{2s}}+\left(1-\alpha \right){\left(\frac{L_D}{L_M}\right)}^2}{\left(1-\alpha \right)+\beta \frac{L_D}{\rho T_{2s}}+\left(1-\alpha \right){\left(\frac{L_D}{L_M}\right)}^2}& \left(4\right)\end{array}$
It should be appreciated that the use of a single transition length L.sub.M in equation (4) causes equation (4) to be effective for uniform pore size samples but does not necessarily address the case of a wide pore size distribution such as in rock formations. It should also be appreciated that in cases involving very large pores and low surface relaxivity where T.sub.2s approaches infinity, equation (4) reduces to

(22) 0$\begin{array}{cc}\frac{D\left(\Delta \right)}{D_0}=\frac{1+{\alpha \left(\frac{L_D}{L_M}\right)}^2}{1+{\left(\frac{L_D}{L_M}\right)}^2}.& \left(5\right)\end{array}$
As 0<α<1 and

(23) $\frac{L_D}{L_M}$
is finite using a fixed L.sub.M, equation (5) does not asymptote to the bulk diffusion coefficient.

(24) In order to permit equation (4) to be effective for a wide pore size distribution and have it approach the bulk diffusion coefficient when T.sub.2s approaches infinity, a parameter Γ is introduced such that L.sub.M is defined by

(25) $\begin{array}{cc}{L}_M=\Gamma *\frac{V}{S}=\Gamma \rho T_{2s}& \left(6\right)\end{array}$
This definition causes L.sub.M to scale with pore size. For large pores, longer diffusion length is required to observe the transition between short-time and long-time diffusion limits. After introducing the new parameter, equation (5) will give the expected bulk diffusion coefficient in very large pores because as T.sub.2s approaches ∞, L.sub.M also approaches infinity. Substituting (6) into (4) now gives the mathematical description of the diffusion coefficient as a function of both diffusion encoding time Δ and relaxation time T.sub.2s with the desired asymptotic behaviors:

(26) $\begin{array}{cc}\frac{D\left(\Delta ,T_{2s}\right)}{D_0}=1-\left(1-\alpha \right)\frac{\beta \frac{L_D}{\rho T_{2s}}+\left(1-\alpha \right){\left(\frac{L_D}{\Gamma \rho T_{2s}}\right)}^2}{\left(1-\alpha \right)+\beta \frac{L_D}{\rho T_{2s}}+\left(1-\alpha \right){\left(\frac{L_D}{\Gamma \rho T_{2s}}\right)}^2}.& \left(7\right)\end{array}$

(27) In one aspect, the generalized model in equation (7) can now be applied to samples with a wide distribution of pore sizes. Each pore size has a distinct diffusion coefficient D(Δ) and T.sub.2s, governed by equations (2) and (3) respectively. For a fixed value of T.sub.2s, equation (7), as a function of diffusion time Δ, reduces to the conventional Padé approximation of time dependent diffusion coefficient whose fitting yields the local surface to volume ratio S/V that is relevant to the NMR measurements. For a specified diffusion time Δ, equation (7) is a function of T.sub.2s. In the special case where T.sub.2s approaches ∞, i.e. very large pores, equation (7) reduces to

(28) $\begin{array}{cc}\frac{D\left(\Delta ,T_{2s}\right)}{D_0}=1-\frac{4}{9\sqrt{\pi }}\frac{1}{\rho T_{2s}}\sqrt{D_0\Delta }& \left(8\right)\end{array}$
which is consistent with the known short-time limit models. See, Zielinski, L., et al., “Restricted Diffusion Effects in Staturation Estimates from 2D Diffusion-Relaxation NMR Maps,” SPE Annual Technical Conference and Exhibition (2010).

(29) Turning now to FIG. 1, an embodiment of a method for determining surface relaxivity of a porous medium is provided. At 10, an NMR tool such as may be used in conjunction with a core analyzer uphole or such as a downhole tool described hereinafter with respect to FIG. 7A or 7B is activated to generate a pulsed field gradient NMR pulse sequence and data (decay signals over time τ) are collected. The pulse sequence is seen in FIG. 2 and may be utilized to obtain a DT.sub.2 distribution containing information about the geometry of a porous medium from which surface relaxivity information may be obtained as described hereinafter. The pulse sequence of FIG. 2 includes on a first channel, two spaced RF pulses of ninety degrees followed by a CPMG sequence of a ninety degree pulse followed by a series of one hundred eighty degree pulses, with signal acquisition prior to and between the one hundred eighty degree pulses, and on a second channel, a gradient pulse between the two spaced RF pulses of ninety degree pulses, a crusher pulse between the second RF pulse and the ninety degree pulse of the CPMG sequence, and another gradient pulse between the ninety degree pulse and the first echo of the CPMG sequence. The separation of the two gradient pulses is the diffusion encoding time Δ. As will be discussed hereinafter, according to one aspect, the pulse sequence of FIG. 2 is repeated for multiple different diffusion times, and for each diffusion encoding time, the pulse sequence is repeated for multiple different encoding strengths (q). In one embodiment, the encoding strengths increase in amplitude from pulse sequence to pulse sequence. In one embodiment, a first pulse sequence (q.sub.0) has gradient pulses having a strength of zero; i.e., they are not present. It will be appreciated that any of many techniques may be used to encode for D(Δ) each with its own advantages and disadvantages with respect to performance and sensitivity to artifacts. See, e.g., Price, W. S., “Pulsed-field gradient nuclear magnetic resonance as a tool for studying translation diffusion: Part 1. Basic theory,” Concepts in Magnetic Resonance, 9(5), pp. 299-336 (1997). Furthermore the pulse sequence components encoding for diffusion and T2 need not be in separate sections in time but may be combined together to simultaneously encode for both. For example, the ordinary CPMG sequence can be made to correlate diffusion and relaxation (T2) in a permanent gradient by, instead of varying a gradient strength, varying the echo time (te) of the entire CPMG. In this way the signal decay rate incorporate both diffusion and T2 information and varying te changes the sensitivity to diffusion and the effective diffusion encoding time (Δ) is the echo time (te). Accordingly, it should be understood that diffusion editing in the context of the present application may encompass (a) diffusion encoding that is a separate sequence from the T2 encoding sequence and (b) diffusion encoding that is at least partially integrated into the T2 encoding sequence.

(30) Returning to FIG. 1, at 20, a DT.sub.2 distribution is determined from the NMR tool signals M(q, τ) resulting from the pulse sequences of different gradient strengths. More particularly, under Gaussian approximation the measured signal by pulsed gradient stimulated echo (PGSTE) followed by CPMG detection is described by

(31) $\begin{array}{cc}M\left(q,\tau \right)=\int \int f\left(T_2,D\right)e^{-{\gamma }^2{g}^2{\delta }^{2}D\Delta }{e}^{-\frac{\tau }{T_2}}dDd{T}_2& \left(9\right)\end{array}$
where f(T.sub.2,D) is the DT.sub.2 distribution, γ is proton gyromagnetic ratio, g is pulsed gradient strength, δ is gradient duration, D is the apparent diffusion coefficient (which includes the effect of restricted diffusion), Δ is diffusion encoding time (of the gradient pulse sequence) and τ is the time in CPMG decay. The encoding strength is characterized by q=γgδ. As previously mentioned, for a given diffusion encoding time, the encoding strength is varied for a series of pulse sequences. This is accomplished by changing the pulsed gradient strength g, and/or by changing the gradient duration δ while keeping Δ fixed. Thus, M(q,τ) are the measured values for discrete values of q and τ and the collected data may be expressed as M(q.sub.i,τ.sub.j).

(32) In some implementations, a 2D fast Laplace inversion algorithm is used to obtain the DT.sub.2 distribution f(T.sub.2, D) (for a particular diffusion encoding time Δ) from the measured values M(q,τ). See, e.g., Venkataraman, L., et al., “Solving Fredhom integrals of the first kind with tensor product structure in 2 and 2.5 dimensions,” Signal Processing, IEEE Transactions Vol. 50 pp. 1017-1026 (2002); and Song, Y.-Q., et al., “T1-T2 Correlation Spectra Obtained Using a Fast Two-Dimensional Laplace Inversion,” Journal of Magnetic Resonance, Vol. 154 pp. 261-268 (2002). More particularly, equation (9) may be rewritten as

(33) $\begin{array}{cc}M\left(q_i,{\tau }_j\right)={\Sigma }_{k,l}f\left(T_{2,k},D_l\right)e^{-q_i^2{D}_l\Delta }{e}^{-\frac{{\tau }_i}{T_{2,k}}}& \left(9a\right)\end{array}$
which represents a linear system of equations that can be solved using a fast numerical Laplace inversion (FLI). The result is a distribution or map of signals as a function indexed by their value of T.sub.2 and D; f.sub.T.sub.2,k.sub.,D.sub.l.

(34) From the DT.sub.2 distribution (also called a “DT.sub.2 map”), according to some implementations, a mean or average diffusion coefficient for a particular Δ is calculated by

(35) $\begin{array}{cc}\stackrel{\_}{D}\left(\Delta ,T_{2s}\right)=\frac{{\Sigma }_{i}f\left(T_{2s},D_i\right)*D_i^2}{{\Sigma }_{i}f\left(T_{2s},D_i\right)*D_i}& \left(10\right)\end{array}$
which results in a diffusion coefficient as function of both T.sub.2s and diffusion encoding time Δ. More particularly, and as seen in FIG. 1 at 30, T.sub.2 indices are converted into T.sub.2s indices for the DT.sub.2 distribution data matrix f.sub.T.sub.2,k.sub.,D.sub.l by relabeling the T.sub.2 index (data axis) as a T.sub.2s index (data axis) utilizing

(36) $\frac{1}{T_{2s}}=\frac{1}{T_2}-\frac{1}{T_{2b}}.$
Then, at 35, the average diffusion coefficient D is found as a function of T.sub.2 according to equation (10). Because magnetization signal decay is usually slowed down at large wave vectors due to the non-Gaussian behavior of restricted diffusion which potentially yields a secondary slow diffusion artifact peak, the D(Δ,T.sub.2s) calculation of equation (10) weights more strongly the large diffusion coefficients. As a result, the effect of non-Gaussian behavior is minimized to give a more accurate determination of the apparent diffusion coefficient. In another embodiment, instead of utilizing equation (10), the average diffusion coefficient is calculated according to

(37) $\begin{array}{cc}\stackrel{\_}{D}\left(\Delta ,T_{2s}\right)=\sqrt{\frac{{\Sigma }_{i}f\left(T_{2s},D_i\right)*D_i^2}{{\Sigma }_{i}f\left(T_{2s},D_i\right)}}.& \left(10a\right)\end{array}$
In some implementations, instead of utilizing equation (10) or equation (10a), the average diffusion coefficient is calculated according to

(38) 0$\begin{array}{cc}\stackrel{\_}{D}\left(\Delta ,T_{2s}\right)=\frac{{\Sigma }_{i}f\left(T_{2s},D_i\right)*D_i}{{\Sigma }_{i}f\left(T_{2s},D_i\right)}.& \left(10b\right)\end{array}$

(39) At 40, the T.sub.2 distribution across diffusion coefficients D.sub.i for a particular diffusion time Δ is calculated according to
g(Δ,T.sub.2s)=Σ.sub.D.sub.if.sub.Δ(T.sub.2s,D.sub.i)  (11)
Then, at 50, a determination is made as to whether data for all diffusion times Δ have been obtained and/or processed. If not, the method is repeated for a new diffusion time by generating a pulsed field gradient NMR pulse sequence with the new diffusion time and collecting data at 10, determining the DT.sub.2 distribution at 20, finding the average diffusion coefficient D at 35 (after relabeling at 30), and calculating the T.sub.2 distribution across diffusion coefficients D.sub.i at 40. It is noted that instead of sequentially generating the pulse sequence and collecting data followed by determining the DT.sub.2 distribution, finding the average diffusion coefficient, etc., some or all pulse sequences may be generated and data collected prior to determining the DT.sub.2 distribution, finding the average diffusion coefficient, etc.

(40) With the T.sub.2 distributions and the average diffusion coefficients D calculated for all Δ, it is then possible to make a determination of the surface relaxivity ρ. More particularly, at 60, a minimization procedure is utilized where equation (7) is fitted to the 2D data set D(Δ,T.sub.2s) in order to extract a determination of the surface relaxivity ρ along with other (unknown) parameters such as the bulk diffusion coefficient D.sub.0, the asymptotic value α=D(∞)/D.sub.0, and parameter Γ=L.sub.m/ρT.sub.2s. In the minimization procedure, the cost function is constructed as follows:
err=Σ.sub.Δ,T.sub.2s{g(Δ,T.sub.2s)*[log D.sub.fit(Δ,T.sub.2s)−log D(Δ,T.sub.2s)]}.sup.2  (12)
where D.sub.fit(Δ,T.sub.2s) is the fitted function of equation (7) whose unknown parameters are chosen to cause the fitted function to most nearly approximate the data set of average diffusion coefficients, thereby causing equation (12) to be a minimum, and g(Δ,T.sub.2s) is the error weighing function, i.e. the T.sub.2 distribution calculated according to equation (11). In the minimization of equation (12), a logarithmic error is used so that it equally weighs the faster and slower diffusion coefficients (due to different pore sizes). In other embodiments, equation (12) is modified to include a non-logarithmic error. In those other embodiments, the slower diffusion coefficient region will have little effect on the fitting.

(41) Alternatively or additionally, in order to generate a function relating Diffusion and T2, from a DT2 map, one can also calculate various weighted means T.sub.2s as a function of diffusion (e.g., T.sub.2s (D.sub.Δ)) to compute the functional relation between D and T.sub.2. For example, the log mean T.sub.2 as a function of D:

(42) $\stackrel{\_}{T_{2s}}\left(D_{\Delta }\right)=10\frac{{\Sigma }_{i}F\left(T_{2s,i},D_{\Delta }\right){\log }_{10}{T}_{2s,i}}{{\Sigma }_{i}F\left(T_{2s,i},D_{\Delta }\right)}$
or the linear mean

(43) $\stackrel{\_}{T_{2s}}\left(D_{\Delta }\right)=\frac{{\Sigma }_{i}F\left(T_{2s,i},D_{\Delta }\right)T_{2s,i}}{{\Sigma }_{i}F\left(T_{2s,i},D_{\Delta }\right)}$

(44) In some implementations, the pulsing at 10, the determining of the DT.sub.2 distribution at 20, the relabeling at 30, the finding of the average diffusion coefficient D as a function of T.sub.2 according to equation (10) at 35, and the calculation of a T.sub.2 distribution at 40 are done for a single encoding time Δ, and equation (7) is fitted to D(T.sub.2s) in order to extract a determination of the surface relaxivity ρ along with other (unknown) parameters such as the bulk diffusion coefficient D.sub.0, the asymptotic value α=D(∞)/D.sub.0, and parameter Γ=L.sub.m/ρT.sub.2s.

(45) According to one aspect, it will be appreciated that the method of FIG. 1 is accomplished utilizing information obtained from NMR data only, and does not require other non-NMR measurements such as measurements of formation parameters such as porosity.

(46) According to one embodiment, multiple determinations of surface relaxivity may be made along a distance or depth in a borehole and the determinations may be displayed as a log according to distance or depth in a borehole. According to another embodiment, the determination of surface relaxivity may be displayed on paper or on a screen.

(47) In order to confirm the method of FIG. 1, four packed glass beads (class IV soda lime spheres) were used as a model system. The glass beads had nominal diameters of 50 μm, 100 μm, 200 μm and 500 μm respectively. The diameter variation was about ±15% and spherical percentage was greater than 85%. The as-received glass beads were cleaned with 1M HCl at room temperature for one hour and then thoroughly washed with deionized water to remove paramagnetic impurities and surface contaminations. The cleaned glass beads were saturated with deionized water whose bulk relaxation time T.sub.2b was 2.3 s. The beads' pack information with respect to bead sizes, porosities and T.sub.2 relaxation times are summarized in Table 1 below.

(48) In accordance with FIG. 1, NMR measurements were carried out on a Magritek Rock Core Analyzer (.sup.1H frequency 2 MHz) at room temperature. The gradient was calibrated using bulk deionized water in 3 mm capillaries to avoid convection. Pulsed gradient stimulated echo (PGSTE) was used for diffusion encoding and the standard CPMG was used for data acquisition as depicted in FIG. 2. The echo time TE between two π pulses was 167 μs. Signal decay due to diffusion in background gradient (<0.02 G/cm) during the CPMG acquisition was negligible.

(49) The measured signal resulting from the NMR pulse sequences was processed according to FIG. 1. Representative fitting results on glass beads are shown in subplot FIGS. 3a-3e (100 μm beads) and FIGS. 4a-4e (500 μm beads). The D(Δ,T.sub.2s) data are displayed in gray scale according to the weighing function g(Δ,T.sub.2s), i.e. T.sub.2 distribution. On each subplot, g(Δ,T.sub.2s) has been normalized for display purposes. Restricted diffusion is clearly seen as diffusion time increases from 40 ms (FIGS. 3a and 4a) to 1 s (FIGS. 3e and 4e), but the T.sub.2 distribution as seen in FIGS. 3f and 4f remains similar for all DT.sub.2 maps except the intensity decreases. Equation (7) is able to fit the diffusion coefficient behavior depending on T.sub.2 relaxation time and diffusion coefficient very well.

(50) The multiple DT.sub.2 global fitting parameters are summarized in Table 1. Bulk liquid diffusion coefficient D.sub.0 varies a little bit because temperature is not strictly controlled and the water diffusion coefficient depends sensitively on temperature. The resultant long time limit diffusion coefficient D(∞)/D.sub.0α falls between 0.59 to 0.66, close to the expected value α=√{square root over (φ)}=0.62 for packed glass beads with porosity φ of about 0.37. It is interesting to note that the parameter Γ stays fairly constant although the bead diameter changed from 50 μm to 500 μm. This result is in agreement with the above suggestion that L.sub.M scales with pore size (and bead size).

(51) TABLE-US-00001 TABLE 1 Sample information and multiple DT.sub.2 global fitting results based on eq. (7). size range porosity D.sub.0 Sample name (μm) φ T.sub.2 (s) ρ (μm/s) (10.sup.−9 m/s.sup.2) α Γ  50 μm bead 45-63 37% 0.27 14.9 2.19 0.66 2.5 100 μm bead  90-125 36% 0.48 16.5 2.18 0.62 2.1 200 μm bead 180-250 37% 0.66 21.8 2.15 0.62 2.1 500 μm bead 425-600 38% 1.19 20.7 2.22 0.59 1.9 Berea Null 23% distribution 19.4 2.15 0.56 2.1 Fontainebleau Null 5.6%  distribution 3.0 2.15 0.62 1.9

(52) To further estimate the reliability of the method, the fitted surface relaxivity ρ was compared to the value calculated by

(53) $\begin{array}{cc}\rho =\frac{d*\varphi }{6*T_{2s}*\left(1-\varphi \right)}& \left(13\right)\end{array}$
where φ is porosity and d is nominal bead size. Here

(54) $\frac{6*\left(1-\varphi \right)}{d*\varphi }$
can be used as a good estimate of surface volume ratio S/V because the bead size is fairly uniform and bead shape is mostly spherical. Since the diffusion length is in the order of micrometers, the surface roughness in nanometer scale should not affect the NMR relevant S/V. The calculated surface relaxivity using the method of FIG. 1 are shown in Table 2 along with single DT.sub.2 fitting results using fixed L.sub.M=100 μm and L.sub.M=20 μm respectively.

(55) TABLE-US-00002 TABLE 2 Comparison between multiple DT.sub.2 global fit and fixed L.sub.M single DT.sub.2 fit Single DT.sub.2 Single DT.sub.2 Nominal Calculated L.sub.M = 100 μm, L.sub.M = 20 μm, bead size ρ (μm/s) Multiple DT.sub.2 Δ = 82 ms, Δ = 82 ms, d by Eq. 13 ρ (μm/s) ρ (μm/s) ρ (μm/s)  50 μm 16.0 14.9 10.1 12.2 100 μm 15.5 16.5 9.8 13.9 200 μm 21.2 21.8 12.6 28.2 500 μm 20.7 20.7 17.1 Null
It is interesting to note that multiple DT.sub.2 fitted surface relaxivity is consistent with the calculated value for all four beads within a ±7% error. However, fixed L.sub.M parameterization on a single DT.sub.2 fitting, while providing a fairly good estimate for a specific bead size, is seen to not satisfy all pore sizes with multiple length scales. For example, while the fitting on 500 μm beads using L.sub.M=100 μm yields ρ=17.1 μm/s, which is only 17% smaller than the calculated 20.7 μm/s using equation (13), on 100 μm beads, it yields ρ=9.8 μm/s which is 37% smaller than the expected surface relaxivity 15.5 μm/s. The misfit is even more serious for a smaller fixed L.sub.M, for instance, L.sub.M=20 μm. For L.sub.M=20 μm, the 500 μm beads could not be fitted with this parameterization because the asymptotic value in equation (5) is significantly below 1.

(56) NMR DT.sub.2 experiments were also carried out on Berea sandstone sample (porosity 23.17%, permeability 1026 mD) and Fontainebleau sandstone sample (porosity 5.63%, permeability 5.15 mD) rock cores where the multiple DT.sub.2 global fitting method of FIG. 1 was applied to determine the surface relaxivities of the samples. FIGS. 5a-5e show the results on the Berea sandstone rock core sample with a wide T.sub.2 distribution indicated in the subplot of FIG. 5f. As diffusion time increases, the smaller T.sub.2 components gradually disappear, shifting the T.sub.2 distribution peak towards larger T.sub.2 values. The fitting of the method yields a surface relaxivity 19.3 μm/s, which is almost twice that of a previously reported value (9.8 μm/s) on a Berea sample with porosity 22% and permeability 380 mD. However, the previously reported value was based on a T.sub.2 distribution peak of 1.0 s while the measured T.sub.2 peak of this experiment was at 0.60 s. Thus, the difference in surface relaxivity might be caused by lithology or pore size differences. In FIGS. 5a-5e, it can be seen that the global fitting based on equation (7) works very well for smaller diffusion time, but an apparent deviation is observed for Δ=500 ms where the D(Δ,T.sub.2s) data lie below the fitted curve for the small and large T.sub.2 region. The D(Δ,T.sub.2s) data has a saddle-like shape, probably due to an internal gradient and the resolution limit of the fast Laplace inversion algorithm on a weak signal. However, this is accounted for by the weighing of the T.sub.2 distribution which focuses more on the well-behaved region. With the method of FIG. 1, the effect of an abnormally small T.sub.2 and a large T.sub.2 region is minimal because the corresponding signal is very weak.

(57) Fontainebleau sandstone results are shown in FIGS. 6a-6e. The T.sub.2 distribution shown in FIG. 6e is between 1 s and 2 s, indicating a relative narrow pore size distribution. However some abnormal features are present on the D(Δ,T.sub.2s) data. For example, on the Δ=40 ms data, the diffusion coefficient decreases with T.sub.2 relaxation time. This could be caused by the resolution of the 2D inversion considering the narrow T.sub.2 distribution and the weak signal due to low porosity. The error due to 2D inversion can be averaged out to give a reasonable fitting result. The fitted surface relaxivity is 3.0 μm/s and other parameters fall into the expected range (Table 1). The small surface relaxivity is in line with the lithology and confirmed by the qualitative observation that its T.sub.2 peak is at 1.36 s, much larger than that of Berea sandstone (0.60 s), even though its pore sizes are smaller than that of Berea. Prior work reported surface relaxivities in the range 5-10 μm/s on a suite of 7 Fontainebleau sandstone samples. Particularly, on a Fontainebleau sample (porosity 6.4%, permeability 5.4 mD) similar to the investigated one, a surface relaxivity 5.1 μm/s was obtained. While the relaxivity determination obtained using the method of FIG. 1 of 3.0 μm/s is smaller than reported value of 5.1 μm/s, the results are consistent with the fact that the reported T.sub.2 distribution peak was at 0.83 s while the T.sub.2 distribution peak in the experiment was at 1.36 s.

(58) It should be appreciated that the method of FIG. 1 may be carried out by any of numerous apparatus. Thus, in one embodiment, the method may be accomplished in an (uphole) laboratory through the use of a core sampler NMR apparatus such as, by way of example only, a Magritek Rock Core Analyzer. In another embodiment, the method may be accomplished utilizing downhole equipment such as shown in FIG. 7A or 7B.

(59) FIG. 7A shows an example nuclear magnetic resonance (“NMR”) wireline well logging instrument 110 disposed in a wellbore 117 drilled through subsurface rock formations 126, 124. The instrument 110 is attached to one end of an armored electrical cable (“wireline”) 118. The cable 118 may be extended into the wellbore 117 and withdrawn therefrom by a spooling device such as a winch 120 of types well known in the art. The cable 118 includes one or more insulated electrical conductors and may include one or more optical fibers to communicate signals between the instrument 110 and a recording unit 122 disposed at the Earth's surface. The recording unit 122 may include a processor or computer having an optional screen or printer type data display, input controls and a data recording device for storage of signals (e.g., NMR measurements) communicated from the well logging instrument 110, as well as for storing or displaying calculated results made from NMR measurements made by the instrument 110 as described hereinafter. The well logging instrument 110 may also include a downhole processor as described hereinafter.

(60) The NMR instrument 110 includes a magnet 112 for inducing a static magnetic field in the formations 124, 126 having a predetermined spatial distribution of magnetic field amplitude. According to one aspect, magnet 112 may be supplemented by an electromagnet configured to impart a selected magnitude gradient field superimposed on the static homogenous field. The NMR instrument 110 may be arranged as disclosed in U.S. Patent Application Publication No. 20140184220, and/or in U.S. Pat. No. 5,796,252 which are both hereby incorporated by reference herein in their entireties. As the instrument 110 is moved along the interior of the wellbore 117, nuclei in the formations surrounding the wellbore are magnetically polarized along the direction of the magnet's 112 field. The instrument 110 also includes an antenna for inducing radio frequency (“RF”) magnetic fields in the formations, and for detecting radio frequency signals induced by NMR phenomena excited in the formations by the static and RF magnetic fields. The particular portion of the formations adjacent to the wellbore from which the NMR signals originate depends on, among other factors, the spatial amplitude distribution of the static magnetic field and the RF frequency used to induce NMR phenomena in the formations. Some magnets may induce a region of substantially homogeneous field amplitude in a particular region in the formations; other types of magnets may induce static fields having a selected amplitude gradient in a particular region of interest. However arranged, the NMR instrument 110 is equipped to generate NMR pulses such as shown in FIG. 2.

(61) Some formations, for example the one illustrated at 124 in FIG. 7A may be permeable and/or contain movable hydrocarbon in the pore spaces thereof. Proximate the wall of the wellbore 117, a portion of the formation 124 may be subjected to sufficient infiltration of the liquid phase of a fluid (“drilling mud”), called “mud filtrate”, used to drill the wellbore 117, that substantially all of the mobile connate fluids in the pore spaces of the formation 124 are displaced by the mud filtrate. Depending on, for example, the fractional volume of pore space (porosity) of the formation 124, and the filtrate characteristics of the drilling mud, the mud filtrate will fully displace all the mobile connate fluids to a depth represented by d.sub.xo in FIG. 7A. The foregoing is referred to as the diameter of the “flushed zone.” Partial displacement of connate fluid is shown extending to a diameter represented by d.sub.i, which is used to represent the diameter of the “invaded zone.” At a certain lateral depth in the formation 124, beyond the diameter of the invaded zone, connate fluid is substantially undisturbed. A quantity of interest in determining possible fluid production in from the formation is the fractional volume of the pore space that is occupied by water (and its complement assumed to be occupied by hydrocarbons). In the uninvaded zone, such fractional volume, called “saturation”, is represented by S.sub.w. Invaded zone and flushed zone water saturations are represented, respectively, by S.sub.i and S.sub.xo.

(62) The example instrument shown in FIG. 7A is only for purposes of explaining a source of measurements that may be used and is not intended to limit the configurations of NMR well logging instrument that may be used to provide measurements of various embodiments hereof. Further, reference to portions of formations that contain hydrocarbon are only for purposes of illustrating general principles of NMR well logging; as certain measurements of NMR properties may be made in formations known to be fully water saturated to simplify calculations of formation properties made from the NMR measurements.

(63) FIG. 7B illustrates a well site system in which an NMR well logging instrument can be conveyed using a drill string or other pipe string for measurement during the drilling of the wellbore, or during other pipe string operations associated with the construction of a wellbore such as circulating, washing, reaming and “tripping.” The well site can be onshore or offshore. In the example system of FIG. 7B, a wellbore 311 is drilled through subsurface formations by rotary drilling in a manner that is well known in the art. Other examples of NMR instruments can be used in connection with directional drilling apparatus and methods. Accordingly, the configuration shown in FIG. 7B is only intended to illustrate a possible source of NMR measurements and is not intended to limit the scope of this disclosure.

(64) A drill string 312 is suspended within the wellbore 311 and includes a bottom hole assembly (“BHA”) 300 proximate the lower end thereof. The BHA 300 includes a drill bit 305 at its lower end. The surface portion of the well site system includes a platform and derrick assembly 310 positioned over the wellbore 311, the assembly 310 including a rotary table 316, kelly 317, hook 318 and rotary swivel 319. The drill string 312 is rotated by the rotary table 316, which is itself operated by well known mechanisms. The rotary table 316 engages the kelly 317 at the upper end of the drill string 312. The drill string 312 is suspended from the hook 318. The hook 318 is attached to a traveling block, through the kelly 317 and the rotary swivel 319 which permits rotation of the drill string 312 relative to the hook 318. As is well known, a top drive system could alternatively be used instead of the kelly 317 and rotary table 316 to rotate the drill string 312 from the surface. The drill string 312 may be assembled from a plurality of segments 325 of pipe and/or collars threadedly joined end to end.

(65) In the present example, the surface system further includes drilling fluid (“mud”) 326 stored in a tank or pit 327 formed at the well site. A pump 329 delivers the drilling fluid 326 to the interior of the drill string 312 via a port in the swivel 319, causing the drilling fluid 326 to flow downwardly through the drill string 312 as indicated by the directional arrow 308. The drilling fluid 326 exits the drill string 312 via water courses, or nozzles (“jets”) in the drill bit 305, and then circulates upwardly through the annulus region between the outside of the drill string and the wall of the borehole, as indicated by the directional arrows 309. In this manner, the drilling fluid 326 lubricates the drill bit 305 and carries formation cuttings up to the surface, whereupon the drilling fluid 326 is cleaned and returned to the pit 327 for recirculation.

(66) The bottom hole assembly 300 of the illustrated example may include a logging-while-drilling LWD) module 320, a measuring-while-drilling (MWD) module 330, a steerable directional drilling system such as a rotary steerable system and/or an hydraulically operated motor such as a steerable motor, and the drill bit 305.

(67) The LWD module 320 is housed in a special type of drill collar, such as, for example, one known in the art, and may contain one or a plurality of known types of well logging instruments. It will also be understood that more than one LWD and/or MWD module may be used, e.g. as represented at 320A. (References, throughout, to a module at the position of LWD module 320 can alternatively mean a module at the position of MWD module 320A as well.) The LWD module 320A typically includes capabilities for measuring, processing, and storing information, as well as for communicating with the surface equipment. In the present embodiment, the LWD module 320 includes an NMR measuring instrument. An example configuration of such instrument is explained above with reference to FIG. 7A.

(68) The MWD module 330 is also housed in a special type of drill collar, such as, for example, one known in the art, and can contain one or more devices for measuring characteristics of the drill string and drill bit. The MWD module 330 further includes an apparatus for generating electrical power for the downhole portion of the well site system. Such apparatus typically includes a turbine generator powered by the flow of the drilling fluid 326, it being understood that other power and/or battery systems may be used while remaining within the scope of the present invention. In the present example, the MWD 330 module may include, for example, one or more of the following types of measuring devices: a weight-on-bit measuring device, a torque measuring device, a vibration measuring device, a shock measuring device, a stick slip measuring device, a direction measuring device, and an inclination measuring device.

(69) The foregoing examples of wireline and drill string conveyance of a well logging instrument are not to be construed as a limitation on the types of conveyance that may be used for the well logging instrument. Any other conveyance known in the art may be used, including without limitation, slickline (solid wire cable), coiled tubing, well tractor, and production tubing.

(70) A recording unit 122A may be disposed at the surface and may include data acquisition, recording, input, control and display devices similar to those of the recording unit shown at 122 in FIG. 7A.

(71) In example methods, measurements of nuclear magnetic resonance (“NMR”) properties of subsurface formations may be made at one or more lateral depths into the formations adjacent to the wellbore. A NMR instrument, as explained above with reference to FIGS. 7A and 7B, can be moved along a wellbore drilled through subsurface formations.

(72) In some examples, NMR measurements may be made using an instrument identified by the trademark MR Scanner which is a trademark of the assignee the present invention. In other examples, the NMR measurements may be made using an instrument identified by the trademark CMR-Plus which is also a trademark of the assignee of the present invention. The NMR instrument, irrespective of type, is generally moved longitudinally along the wellbore and a record with respect to depth in the wellbore is made of the NMR properties of the various formations. The foregoing identified MR Scanner instrument, in particular, can make measurements of NMR properties of the formations at a plurality of different, defined lateral depths of investigation. The lateral depths of investigation for the foregoing instrument are about 1.5 inches (3.8 cm), 2.7 inches (6.9 cm) and 4 inches (10.2 cm) from the wall of the wellbore. As explained above, the lateral depth of investigation of any particular NMR measurement is defined by the spatial distribution of the amplitude of the static magnetic field and the frequency of the RF magnetic field used to excite NMR phenomena. The example instruments described herein are not limitations but are provided as illustrative examples.

(73) According to some examples, the processor or computer recording system 122 that is provided may include one or more processors located uphole and/or downhole and can be used to implement or perform any of the methods and processes for analyzing described hereinafter. Where, for example, the method is to be accomplished in a laboratory through the use of a core sampler NMR apparatus, it will be appreciated that any processor or computer recording system may be located uphole. The terms “processor” and “computer recording system” (hereinafter “processing system”) should not be construed to limit the embodiments disclosed herein to any particular device type or system. The processing system may be, for example, one or more laptop computer, one or more desktop computer, and/or one or more mainframe computer. The processing system may also include a processor (e.g., a microprocessor, microcontroller, digital signal processor, or general purpose computer) for executing any of the methods and processes hereinafter described. In this regard, the term “processor” as used herein should be understood, unless indicated to the contrary, to encompass single or multiple processing units disposed in any suitable location or locations. The processing system may further include a memory such as a semiconductor memory device (e.g., a RAM, ROM, PROM, EEPROM, and/or Flash-Programmable RAM), a magnetic memory device (e.g., a diskette and/or fixed disk), an optical memory device (e.g., a CD-ROM), a PC card (e.g., PCMCIA card), or other memory device. This memory may be used to store, for example, fundamental equation and/or instructions for performing the processes described above.

(74) Any of the methods and processes hereinafter described can be implemented as computer program logic for use with the processing system. The computer program logic may be embodied in various forms, including a source code form or a computer executable form. Source code may include a series of computer program instructions in a variety of programming languages (e.g., an object code, an assembly language, or a high-level language such as C, C++, or JAVA). Such computer instructions can be stored in a non-transitory computer readable medium (e.g., memory) and executed by the processing system. The computer instructions may be distributed in any form as a removable storage medium with accompanying printed or electronic documentation (e.g., shrink wrapped software), preloaded with a computer system (e.g., on system ROM or fixed disk), or distributed from a server or electronic bulletin board over a communication system (e.g., the Internet or World Wide Web).

(75) As used herein (including in the claims), a recitation in the general form of “at least one of [A] and [B]” should be construed as disjunctive. That is, this recitation should be construed to include each of the following: (i) both [A] and [B] are present; (ii) [A] is present, but [B] is not, and (iii) [B] is present, but [A] is not.

(76) There have been described and illustrated herein several embodiments of methods of determining surface relaxivity of a porous medium. While particular embodiments and aspects have been described, it is not intended that the disclosure be limited thereto, and it is intended that the claims be as broad in scope as the art will allow and that the specification be read likewise. Thus, while certain equations for obtaining an average diffusion coefficient were described, other equations could be utilized. Similarly, while a particular minimization equation was described, other minimization equations or techniques could be utilized. It will therefore be appreciated by those skilled in the art that yet other modifications could be made. Accordingly, all such modifications are intended to be included within the scope of this disclosure as defined in the following claims. In the claims, means-plus-function clauses, if any, are intended to cover the structures described herein as performing the recited function and not only structural equivalents, but also equivalent structures. It is the express intention of the applicant not to invoke 35 U.S.C. §112, paragraph 6 for any limitations of any of the claims herein, except for those in which the claim expressly uses the words ‘means for’ together with an associated function. Moreover, it should be understood that any of the features described herein may be provided in any suitable combination.