Method of estimating permeability using NMR diffusion measurements

Abstract

This invention is useful for determining the permeability of a geological formation using .sup.1H NMR diffusion measurements acquired in the laboratory and using downhole .sup.1H NMR well logging. The current technology for obtaining formation permeability downhole using NMR is not adequate for low-permeability, unconventional source rock formations with high organic content. This new method uses laboratory .sup.1H NMR diffusion measurements for creating continuous downhole well logs of the mobile-hydrocarbon permeability of the hydrocarbon-filled pore space of downhole geological formations.

Claims

1. A method of mapping, in at least a vertical direction, subsurface permeability as a function of downhole-location, the method comprising: a. acquiring a set CORE_SET of N core samples {CS.sub.1, CS.sub.2 . . . CS.sub.N} where CS.sub.i is the i.sup.th core sample, and N is a positive integer, and 1?i?N, wherein each core-sample CS.sub.i is extracted from a respective downhole extraction-location EXTRACTION-LOC(CS.sub.i) within the subsurface; b. for each core sample CS.sub.i of CORE_SET, respectively performing the following set of measurements: (i) a respective lab-NMR measurement of a respective porosity ?.sub.M(CS.sub.i) describing the porosity of movable hydrocarbons; (ii) a respective lab-NMR measurement of one or more pore-dimension properties; and (iii) a respective lab measurement of the permeability k(CS.sub.i) of the core sample CS.sub.i; c. for each core sample CS.sub.i of CORE_SET, estimating a respective downhole hydraulic tortuosity ?.sub.hy(EXTRACTION-LOC(CS.sub.i)) representing the downhole tortuosity at the extraction-location from which core sample CS.sub.i was extracted; d. computing, from the results of steps (b) and (c), an at least 1-D BTR-map of subsurface BTR (body throat ratio) as a function of subsurface location; e. for a plurality of downhole locations, downhole-NMR-logging the following parameters: (i) subsurface porosity; and (ii) a subsurface pore-size parameter; f. computing an at least 1-D map of subsurface permeability as a function of downhole-location from all of the following: (i) the at least 1-D BTR-map; (ii) the results of the downhole-NMR-logging; and (iii) the results of step (c).

2. The method of claim 1 wherein the lab-NMR-measured pore-dimension property of step b(ii) is a surface-to-volume ratio.

3. The method of claim 1 wherein the lab-NMR-measured pore-dimension property of step b(ii) is a pore diameter.

4. The method of claim 1 wherein step (d) comprises computing, from the results of steps (b) and (c), a respective laboratory BTR(CS.sub.i) value for each core sample CS.sub.i.

5. The method of claim 4 wherein the respective laboratory BTR(CS.sub.i) value is computed using a modified Carman-Kozeny relationship.

6. The method of claim 1, further comprising drilling one or more horizontal or vertical wells in location(s) determined by the content of the at least 1-D map of subsurface permeability.

7. The method of claim 1 further performing at least one of the following in location(s) determined by the content of the at least 1-D map of subsurface permeability: (i) casing and perforating; and (ii) deploying a pump in the well.

8. The method of claim 1 wherein the core samples of CORE_SET are obtained from a selection of depths covering different lithologies (i.e. rock types).

9. The method of claim 1 wherein step (b)(i) comprises saturating, in the laboratory, each core sample with decane and step b(ii) comprises saturating, in the laboratory, each core sample with methane.

10. The method of claim 1 wherein step (b)(i) comprises measuring, in the laboratory, D-T.sub.2, on core samples of CORE_SET to determine a respective movable hydrocarbon fluid porosity ?.sub.M(CS.sub.i) i.e. porosity above T.sub.2>T.sub.2cutoff.

11. The method of claim 1 wherein step (b)(ii) comprises measuring, in the laboratory, restricted diffusion (D/D.sub.0), where D.sub.0 is bulk diffusion versus diffusion length $(L_D=\sqrt{D_0{t}_{?}},$ where t.sub.? is the diffusion encoding time.

12. The method of claim 11 wherein step b(ii) comprises fitting D/D.sub.0 versus L.sub.D using Pad? fit to determine pore body size (d), heterogeneity length scale (L.sub.M), and tortuosity (?) for each core.

13. The method of claim 1 wherein step b(iii) is performed by setting each core sample CS.sub.i in a Hassler sleeve and applying confining stress on the sample CS.sub.i and then flowing fluid (either gas or liquid) through each sample CS.sub.i to produce a pressure drop across the sample CS.sub.i and measuring the fluid flow rate.

14. The method of claim 1 wherein step (c) comprises subjecting each core sample CS.sub.i to laboratory NMR.

15. The method of claim 1 wherein step (c) is performed by at least one of (i) downhole NMR and (iii) downhole electrical resistivity measurements, for example, at EXTRACTION-LOC(CS.sub.i).

16. The method of claim 1 wherein step (f) comprises employing at least one of the modified Carman-Kozeny equations.

17. A method comprising performing NMR restricted diffusion measurements on a core with connate water present by: a) pressure-saturating the core with two hydrocarbons (sequentially high-pressure methane and decane, or high-pressure methane and high-pressure butane), b) measuring NMR restricted diffusion of the two hydrocarbons under pressure in the core to obtain their restricted diffusivity in the hydrocarbon-bearing porosity of the core, c) estimating the tortuosity and surface-to-volume ratio of the hydrocarbon-filled pore space of the core, and d) using the tortuosity, surface-to-volume ratio, and porosity occupied by the hydrocarbon fluids, and estimating the pore body-to-pore-throat ratio for the lithology, determining the permeability of the mobile hydrocarbons in the core.

18. The method of claim 17 wherein the tortuosity and pore size are determined by applying a Pad? fit to restricted diffusivity NMR measurements of the two hydrocarbons.

19. The method of claim 17 wherein the two hydrocarbons comprise methane and at least one other hydrocarbon such as decane or butane.

20. The method of claim 17 wherein the equation for the permeability of the mobile hydrocarbons in the core is given by the equation $k=\frac{?d^2}{32?{\mathrm{BTR}}^2}.$

21. A method comprising obtaining a downhole wireline permeability log by measuring NMR diffusion-T.sub.2 (D-T.sub.2), and permeability on a selection of cores of different lithologies (i.e. rock types), then using the core-calibrated parameters with D-T.sub.2 well logs to obtain a permeability well log, the method comprising the following steps: a. obtaining cores from a selection of depths covering different lithologies, b. measuring D-T.sub.2 on cores to determine movable fluid porosity (?.sub.M), i.e. porosity above T.sub.2>T.sub.2cutoff, and restricted diffusion (D/D.sub.0), where D.sub.0 is bulk diffusion versus diffusion length $(L_D=\sqrt{D_0{t}_{?}},$ where t.sub.? is the diffusion encoding time) using a light hydrocarbon and a heavier hydrocarbon to cover a range of L.sub.D, c. fitting D/D.sub.0 versus L.sub.D using Pad? fit to determine pore body size (d), heterogeneity length scale (L.sub.M), and tortuosity (?) for each core, d. determining the cementation exponent (m) for each lithology using $?={?}_M^{1-m},$ e. determining the heterogeneity length-scale (L.sub.M) for each lithology, f. measuring core permeability (k) and determining the body-throat ratio BTR for each lithology from the Carman-Kozeny equation $k=\frac{{?}_M{d}^2}{32?BT{R}^2},$ g. measuring D-T.sub.2 log to determine movable fluid porosity (?.sub.M, log), i.e. porosity above T.sub.2>T.sub.2cutoff, and D.sub.log/D.sub.0 versus T.sub.2, h. determining tortuosity log (?.sub.log) using ${?}_{log}={?}_{M,\log }^{1-m},$ where m is known for each lithology from step (d), i. fitting D.sub.log/D.sub.0 versus T.sub.2 using Pad? fit to determine pore-body size (d.sub.log), with ?.sub.log known from step (h), and L.sub.M fixed for each lithology from step (e), and j. determining downhole permeability (k.sub.log) from $k_{log}=\frac{{?}_{M,\log }{d}_{\log }^2}{32{?}_{\log }BT{R}^2},$ with BTR known for each lithology from step (f).

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1: Diffusion-T.sub.2 (a.k.a. D-T.sub.2) unipolar stimulated-echo pulse sequence (Tanner 1970, Mitchell 2014) used to measure restricted diffusion as a function of diffusion time t.sub.?. Trapezoidal gradient encoding pulses G (green) are shown with duration t.sub.? (time from half-height to half-height) and ramp-time ?, along with 90? (thin) and 180? (thick) RF pulses (blue) of phase ?.sub.1,2,3,4,a, and CPMG echo trains with echo-spacing t.sub.E. Not shown is a set of 6 woodpecker gradient pulses to stabilize the Eddy currents before the first RF pulse.

(2) FIG. 2: A schematic of the NMR spectrometer with laboratory apparatus for fluid saturation used in this invention.

(3) FIG. 3: The bulk diffusivity D vs. viscosity/temperature ?/T of C1 (methane).

(4) FIG. 4: Illustration of a random walk of a light hydrocarbon molecule in the core sample used in this study, in the short-time limit and the tortuosity limit, at 913 m and 920 m.

(5) FIG. 5: T.sub.2 distribution of the as-received (ASREC), decane saturated (C10), methane-saturated deuterated core C1(D.sub.2O), and deuterated core (D.sub.2O), top is (913 m) and bottom is (920 m). pu1 means assuming the HI=1.

(6) FIG. 6: D-T.sub.2 correlation map of the C4(D.sub.2O) and C1(D.sub.2O) saturated core at 913 m.

(7) FIG. 7: D-T.sub.2 correlation map of the C10 and C1 saturated core at 920 m.

(8) FIG. 8: Diffusion projections from D-T.sub.2 correlation map of the C10 and C1 saturated core at 920 m, normalized to bulk diffusion, with the 2D peak values shown as the data points.

(9) FIG. 9: Diffusion length (L.sub.D) against normalized restricted diffusivity (D/D.sub.0). The dots are the points from the 2D peak of the D-T.sub.2 maps in the region C. Top: C4(D.sub.2O) and C1(D.sub.2O) are shown for 913 m. Bottom: C10 and C1(D.sub.2O) are shown for 920 m. The solid black line is the best fit using the Pad? equation. The dashed horizontal line shows the tortuosity limit. To do the parameter estimation, non-linear curve-fitting (Isqcurvefit) in MATLAB is applied. The loss function is the minimum square error on a log scale.

(10) FIG. 10: Pad? fit in Eq. (14) plotted on the D-T.sub.2 map for C10 where L.sub.D=5 ?m, T.sub.2B=1.33 s, L.sub.M=6.1 ?m, ?=149, and ?=0.92 ?m/s. Also shown is the (poor) result using the electrical tortuosity instead ?.sub.E=5.5, and ?=0.46 ?m/s. Table 3: Fixed parameter settings in the error analysis.

(11) FIG. 11: simulated C1 and C10 data from the Pad? equation. The solid curve is the Pad? equation using Table 3 parameters. It shows data points of C1 and C10 mean value used in the Pad? equation.

(12) FIG. 12: Diffusion length (L.sub.D) against normalized restricted diffusivity (D/D.sub.0) for C10 and C1(D.sub.2O) are shown for 920 m, with additional data compared to FIG. 9

(13) FIG. 13: Cumulative (top) and differential (bottom) distribution in pore-throat diameter from MICP at 913 m and 920 m.

(14) FIG. 14: BTR for 913 m and 920 m, using the entry and mode of the MICP pore-throat distributions. Table 4: Body-throat ratios for different rock types according to (Muller-Huber 2016) *, and the ATR Rock Catalog **.

(15) FIG. 15: Permeability estimates for two organic rich chalks from 913 m and 920 m depths and two Austin chalk samples compared to permeability measurements.

(16) FIG. 16: Illustration of flow paths comprised of a sequence of cone frustrums of pore-body size d and pore-throat size d.sub.throat=d/BTR.

(17) FIG. 17: The value of M as a function of BTR for the permeability equation.

(18) FIG. 18: Diffusion length (L.sub.D) against normalized restricted diffusivity (D/D.sub.0) for 100% H.sub.2O and 100% C1 saturated Indiana limestone.

(19) FIG. 19: Example of decay in magnetization M(b) due to diffusion encoding parameter

(20) $b=q^2\left(t_{?}-t_{?}/3+?\right)$
for an Austin chalk core at S.sub.wirr with C10(D.sub.2O) and L.sub.D=4.9 ?m. Filled symbols show the fitting region.

(21) FIG. 20: Example of restricted diffusion D/D.sub.0 versus L.sub.D determined by the 2D peak (i.e. mode) of the D-T.sub.2 map versus the initial slope, and the resulting Pad? fit, for an Austin chalk core at S.sub.wirr with C10(D.sub.2O) or C1(D.sub.2O).

(22) FIG. 21: D-T.sub.2 correlation map of brine-saturated (SW1) Indiana Limestone.

(23) FIG. 22: Example of restricted diffusion D/D.sub.0 versus L.sub.D determined by the 2D peak (i.e. mode) of the D-T.sub.2 map or the initial slope, and the resulting average tortuosity determined from C1(D.sub.2O), for an Austin chalk core at S.sub.wirr with C10(D.sub.2O) or C1(D.sub.2O).

(24) FIG. 23: Hydrogen Index (HI) comparison between water (unity), CH.sub.4, H.sub.2, and HD, all at 30? C. We assume HD and H.sub.2 have the same molar density at a given pressure and temperature. The hydrogen indexes of fluids other than water are calculated by comparing the proton molar densities with that of water at the same temperature and pressure.

(25) FIG. 24: Example of decay in magnetization M(b) due to diffusion encoding parameter

(26) $b=q^2\left(t_{?}-t_{?}/3+?\right)$
where q=?t.sub.?G, for an Austin chalk core at S.sub.wirr with C1 using the default pulse-parameters in the table.

(27) FIG. 25: Example of decay in magnetization M(b) due to diffusion encoding parameter

(28) $b=q^2\left(t_{?}-t_{?}/3+?\right)$
where

(29) 0$q=?t_{?}G$
for an Austin chalk core at S.sub.wirr with C1 using the optimized pulse-parameters in the table.

(30) FIG. 26: Comparison of diffusion projection D using default versus optimal pulse sequence parameters.

(31) FIG. 27: Comparison of D-T.sub.2 maps using default versus optimal pulse sequence parameters.

DETAILED DESCRIPTION OF EMBODIMENTS

(32) Pore-Size and Tortuosity from NMR Diffusion-T.sub.2 Measurements

(33) In bulk fluid, the molecule diffuses with its bulk diffusivity (D.sub.0). However, with the pore wall's restriction, the measured diffusivity (D) for the pore fluid will be smaller than the bulk diffusivity. As a result, the normalized restricted diffusivity (D/D.sub.0) contains information about the geometry of the connected pore space, including tortuosity (?) and pore size (d), which along with porosity (?) are useful for estimating permeability.

(34) In (Chen 2019, Chen 2019b, Wang 2020) we reported on the diffusive tortuosity and pore size of a low permeability unconventional organic-rich chalk formation at a depth of 920 m. NMR D-T.sub.2 correlation maps of C10 and C1 saturated cores in the presence of connate water are used to determine d and ? in the hydrocarbon-bearing pore space, with connate water present in the micritic pores.

(35) The Pad? fit to the restricted diffusion (D/DO) data as a function of L.sub.D (measured over a large range from 5 ?m to 166 ?m) yields a tortuosity of ?=149 and a pore-body size of d=4.9 ?m for the hydrocarbon-bearing pore space at a depth of 920 m. Error analysis shows ?4.3% accuracy for d and ?1.3% accuracy for ? based on 1% synthetic noise in the simulated measured diffusivity and diffusion-length data. Adding additional data from different hydrocarbons (C2, C3, C4, C5, each at a single L.sub.D) indicates a ?10% change in d and a ?1% change in ?.

(36) The same Pad? fit as a function of T.sub.2 agrees well with the D-T.sub.2 map for C10, with surface relaxivity ?=0.92 ?m/s as the only additional parameter at a depth of 920 m. This agreement helps validate our technique, and most importantly allows for the invention to be used on downhole logs where one cannot change the hydrocarbon-bearing pore space with different fluids.

(37) Experimental

(38) The equation (Tanner, 1970; Mitchell et al., 2014) which describes the magnetization decay after applying the Unipolar Stimulated-Echo pulse sequence is:

(39) $M\left(q,{\mathrm{nt}}_E\right)=\underset{D,T_2}{.Math.}{e}^{-n{t}_E/T_2}{e}^{-q^2\left(t_{?}-t_{?}\right)/3+?)D}P\left(D,T_2\right){?}_D{?}_{T_2}$

(40) where q=?Gt.sub.? is the wavevector for diffusion encoding, ?/2?=42.57 MHz/T is the gyro-magnetic ratio for .sup.1H, ?=?.sup.3/30t.sub.?.sup.2??.sup.2/6t.sub.? is the (small) effect of the ramp time for the gradient pulses, and n is the echo number. P(D, T.sub.2) is the probability density function (i.e. D-T.sub.2 map) normalized as such ?.sub.D,T.sub.2 P(D, T.sub.2) ?.sub.D?T.sub.2=46D.sub.T.sub.2, where ?.sub.DT.sub.2 is the total porosity from D-T.sub.2, and the constants ?.sub.D=log.sub.10(D.sub.i+1)?log.sub.10(D.sub.i), and ?.sub.T.sub.2=log.sub.10(T.sub.2,i+1)?log.sub.10(T.sub.2,i) are the logarithmic bin spacings of the 2D map. The porosity ?.sub.DT.sub.2 is typically less than the total porosity ?.sub.T.sub.1.sub.T.sub.2 from T.sub.1-T.sub.2 maps due to T.sub.1 and T.sub.2 decay during the dead time t.sub.D (see FIG. 1) as such:

(41) ${?}_{D{T}_2}=\underset{T_1,T_2}{.Math.}{e}^{-2{?}_{se}/T_2}{e}^{-\left(t_{?}-{?}_{se}\right)/T_1}p\left(T_1,T_2\right){?}_{T_1}{?}_{T_2}$

(42) where P(T.sub.1, T.sub.2) is the probability density function i.e. the T.sub.1-T.sub.2 normalized as ?.sub.T.sub.1.sub.,T.sub.2 P(T.sub.1, T.sub.2) ?.sub.T.sub.1?.sub.T.sub.2=?.sub.T.sub.1.sub.T.sub.2. Partial compensation of the loss in the signal from D-T.sub.2 can be recovered (numerically) by boosting the P(D, T.sub.2) map as such:

(43) $P_B\left(D,T_2\right)=P\left(D,T_2\right)B\left(T_2\right)$

(44) where the boosting factor B(T.sub.2) at each T.sub.2 bin is:

(45) $B\left(T_2\right)=\min \left\{e^{2{?}_{se}/T_2}{e}^{\left(t_{?}-{?}_{se}\right)/T_{1LM}},3\right\}$

(46) T.sub.1LM(T.sub.2) is defined as the log-mean T.sub.1LM at each T.sub.2 bin, which can be computed from P(T.sub.1, T.sub.2). In cases where P(T.sub.1, T.sub.2) data is not available, one can use a fixed T.sub.1=1.5 T.sub.2, i.e. a constant T.sub.1/T.sub.2=1.5, for instance. The amount of boosting to B(T.sub.2)?3 is limited to avoid boosting noise in the P(D, T.sub.2) map below T.sub.2<T.sub.2,Dcut. In the case of a dead time of t.sub.D=25 ms, a boosting limit of 3 corresponds to T.sub.2,Dcut?7 ms, where P(D, T.sub.2) is mostly noise below T.sub.2<T.sub.2,Dcut.

(47) TABLE-US-00001 TABLE 1 ? t.sub.? ?.sub.se t.sub.? t.sub.D G t.sub.E (ms) (ms) (ms) (ms) (ms) (G/cm) (ms) 0.2 9.0 10 15 110 25 120 0 43 0.2 ?.sub.1 0202020213131313 ?.sub.2 0022002200220022 ?.sub.3 0000222200002222 ?.sub.4 1111111100000000 ?.sub.5 2002022013313113 Top: D-T.sub.2 pulse parameters defined in FIG. 1. Not listed is the spoiler gradient pulse (during the T.sub.1 decay portion) of duration 1 ms and strength 1 G/cm. Bottom: RF phases ?.sub.1,2,3,4,? defined in multiples of ?/2 (i.e. a value of 3 corresponds to a phase of 3?/2).

(48) Note that the boost B(T.sub.2) does not affect the mode of the P(D, T.sub.2) map in region C, as such the results to the Pad? fit shown below are not affected.

(49) The unipolar stimulated-echo pulse sequence developed by Tanner 1970, shown in FIG. 1, has 7 intervals of diffusion encoding, with default parameters listed in Table 1. Another pulse sequence often used in porous media is the bipolar stimulated-echo sequence developed by Cotts et al. 1989, which has 13 intervals of diffusion encoding. While the Cotts 13 sequence is more complex to implement than the Tanner 7 sequence, the Cotts 13 sequence greatly reduces the effects of internal magnetic field gradients in porous media. Internal gradients can distort the diffusion measurement, especially at high magnetic fields, and especially for porous media with large amounts of paramagnetic impurities, such as sandstones. Other variations exist such as the Cotts 17 sequence etc., which can also be used.

(50) A Geospec2 rock-core analyzer by Oxford Instruments was used to make the NMR measurements. A frequency of 2.3 MHz for .sup.1H was used, which is similar to the NMR downhole logging tools. An Oxford NMR overburden cell was used for high-pressure saturation measurements. The NMR spectrometer is equipped with both a radiofrequency probe as well as a gradient coil for making diffusion measurements.

(51) 1D T.sub.2 were acquired using CPMG (Carr-Purcell-Meiboom-Gill) pulse sequences with echo spacing of t.sub.E=0.2 ms. NMR 2D D-T.sub.2 data were acquired using the unipolar stimulated-echo sequence shown in FIG. as a function of diffusion evolution time t.sub.?. The inversion algorithm (Venkataramanan et al., 2002) is used to generate T.sub.2 distributions and D-T.sub.2 maps.

(52) Twin core plugs were taken from a well in the Golan Heights in northern Israel at a depth of 920 m, corresponding to the late Cretaceous upper Ghareb formation. Detailed information about the core samples can be found in the (Chen et al., 2019), including the core plug and well information. The core-plug diameters are 25 mm and the lengths are 48 mm, which are suitable for the NMR overburden cell.

(53) An ISCO pump shown in FIG. 2 was used to inject the hydrocarbon, and a hand pump was used to maintain a net overburden pressure of 1000 psi. The pore pressure was 1200 psia and the spectrometer temperature 30? C. Heat shrink (PTFE) tubing was used to separate the core and the confining fluid. The end piece of the NMR core holder was made of PEEK (polyether ether ketone), which has minimal NMR signal at t.sub.E=0.2 ms.

(54) The bulk properties of hydrocarbon are shown in Table 2. Viscosity and density are obtained from the REFPROP NIST database. The HI (NMR Hydrogen Index) was calculated based on the density and the amount of hydrogen in each molecule. Diffusivity of C1 is from Chen et al., 2019 shown in FIG. 3.

(55) TABLE-US-00002 TABLE 2 Bulk properties of water and hydrocarbons under the experimental conditions. Experiment conditions: 30? C., 1200 psia Viscosity Diffusivity Fluid HI (cP) (?m.sup.2/ms) Water (H.sub.2O) 1.00 0.797 2.3 Methane (C1) 0.13 0.013 250 Decane (C10) 1.02 0.869 1.6

(56) The as-received core contains connate water. As such the C10 saturated core contained both C10 and connate water. A twin core was used for C1 saturation, which was pretreated with D.sub.2O brine (8000 ppm NaCl D.sub.2O brine solution) before C1 saturation to avoid interference from water signal; this core was labeled C1(D.sub.2O).

(57) T.sub.2 Distributions

(58) The rock matrix mainly consists of calcite and kerogen, which are not detected by our NMR measurements. The detectable T.sub.2 range is separated into three regions, A, B and C as shown in FIG. 4 and FIG. 5. FIG. 4 is the illustration of the core sample. Region A (0.2 ms?3 ms) is mainly comprised of bitumen. During the hydrocarbon saturation, the bitumen's viscosity decreases due to hydrocarbon dissolution, which will make the bitumen relax slower, and become more detectable (Singer et al., 2016; Chen et al., 2017). Region B (3 ms?37 ms), which is the micropore region, is mainly comprised of the micritic calcite. Region B contains the connate water in the reservoir with unknown composition. Region C (>37 ms) is comprised of light hydrocarbon and some connate water. Region B is water wet and region C is mix-wet (Chen et al., 2019). During the saturation, some water in the region C is replaced by hydrocarbon. When the diffusion time is short (L.sub.D<d.sub.p) the diffusion is less restricted. When the diffusion time becomes longer (L.sub.D>>d.sub.p), molecules travel to more areas. When the time becomes long enough, the restricted diffusion will approach the limit D.sub.?.

(59) The restricted diffusivity of C10 lies in the short diffusion-length regime (L.sub.D<d.sub.p), which determines the surface-to-volume ratio S/V and pore-body size d.sub.p. The restricted diffusivity of high-pressure C1 lies in the long diffusion-length regime (L.sub.D>>d.sub.p), which determines the tortuosity ?.

(60) For each saturation state, C10 or C1(D.sub.2O), multiple t.sub.? is used in the pulse sequence in order to measure diffusivity at multiple values of L.sub.D, thereby improving the accuracy of pore-body size d.sub.p and tortuosity ? from the Pad? fit.

(61) FIG. 5 shows the T.sub.2 distributions at 913 m (core contains less bitumen) and at 920 m (core contains more bitumen). Region C is mainly comprised of hydrocarbon-filled space, which is the region of interest for mobile hydrocarbons. Region C is mainly hydrocarbon saturated, with some connate water. At higher hydrocarbon viscosity (C10 or C4), we expect some of the connate water to be driven out in region C during saturation. We choose region C as the representative NMR porosity in region C since there is little or no diffusion coupling with the connate water in region B. The porosity in region C is found to be 14.6 pu at 913 m, and 13.0 pu at 920 m, i.e. they are comparable. The slightly larger region C porosity at 913 m may be due to there being less bitumen present than at 920 m (see FIG. 5).

(62) Tortuosity is impacted by the hydrocarbon-filled porosity in region C. Hence to properly study the tortuosity of region C, at 920 m the C1 and C10's saturation in region C should be the same. After adjusting for HI=0.13 of C1, the C1 porosity in region C is ?12 pu which is similar to C10 (13.0 pu), and therefore gives confidence in our value of tortuosity.

(63) Region B mostly contains connate water. At 920 m the peak T.sub.2 value (i.e. mode) for C10 region B is shorter than the as-received. One hypothesis is that C10 replaces some connate water in region B, thereby reducing the connate water in region B and increasing S/V for connate water. Another observation is that NMR porosity in region B for C10 at 920 m is less than as received ASREC. This is because the C10 that entered the micritic pore has slower surface relaxation than brine, and therefore its NMR signal appears in region C.

(64) The shorter T.sub.2 values in region C for 913 m are due to smaller pores and/or larger surface relaxivity than for 920 m.

(65) D-T.sub.2 Correlation Map

(66) D-T.sub.2 measurements (FIG. 6 and FIG. 7) were optimized to get the restricted diffusion information of the hydrocarbon in the region C. Diffusion-T.sub.2 (a.k.a. D-T.sub.2) unipolar stimulated-echo pulse sequence is used to encode the diffusion information into the magnetization decay. After applying a 2D inversion (Venkataramanan et al., 2002), D-T.sub.2 map is acquired to correlate diffusion and T.sub.2.

(67) FIG. 6 shows two representative D-T.sub.2 maps for C4(D2O) and C1(D2O) at 913 m. The y-axis is the restricted diffusivity and the x-axis is the T.sub.2. The porosity is displayed perpendicular to the page, which is shown as contours. The dashed horizontal line represents the bulk diffusivity, which for C1 (250 ?m.sup.2/ms) is out of range of the plot. The oil correlation time (Lo et al., 2002), is shown in the diagonal dash line. Each plot overlays D-T.sub.2 measurements for the same saturated core but with different diffusion lengths listed in the legend. The subplots on the top and right are the T.sub.2 and diffusion projections, respectively. In the top subplot, dashed lines are 1D T.sub.2 distributions from CPMG, and solid lines are from D-T.sub.2 projection. Porosity shown in the legend is the total porosity from 1D T.sub.2 measurement and pu1 means assuming the HI=1.

(68) The restricted diffusion phenomenon is observed in both C4 and C1. Restricted diffusion will change with diffusion length for each hydrocarbon. Different diffusion lengths L.sub.D are acquired by changing the diffusion evolution time t.sub.? in the pulse parameter setting. D-T.sub.2 measurements lose all the region A's signal, and lose most of the region B's signal, due to finite diffusion encoding dead time. The T.sub.2 distribution in region C from the D-T.sub.2 map matches the 1D T.sub.2 measurement in region C (>37 ms). To retrieve the restricted diffusion information of the hydrocarbons, the 2D peak value (i.e. the mode) from the D-T.sub.2 maps is selected. FIG. 7 shows C10 and C1(D2O) at 920 m (more bitumen), which shows similar trends. For completeness, FIG. 8 shows the diffusion projections from FIG. 7, normalized to bulk diffusion, with the 2D peak values shown as the data points.

(69) Pad? Fit

(70) The restricted diffusion is a function of bulk diffusivity (D.sub.0) of the molecule, diffusion time (?.sub.?), surface-to-volume ratio (S/V) of the pore, heterogeneity length-scale (L.sub.M), and diffusive tortuosity (?) (Latour et al., 1993). The Pad? equation describes the relation between these variables (Hurlimann et al., 1994) as such:

(71) $\frac{D}{D_0}=1-\left(1-\frac{1}{?}\right)\frac{\frac{4}{9\sqrt{?}}\frac{S}{V}{L}_D+\left(1-\frac{1}{?}\right)\frac{L_D^2}{L_M^2}}{\frac{4}{9\sqrt{?}}\frac{S}{V}{L}_D+\left(1-\frac{1}{?}\right)\frac{L_D^2}{L_M^2}+\left(1-\frac{1}{?}\right)}$

(72) Where the diffusion length is defined as:

(73) $L_D=\sqrt{D_0{t}_{?}}$

(74) If assuming cylindrical geometry of diameter d:

(75) $\frac{S}{V}=\frac{2}{r}=\frac{4}{d}$

(76) the Pad? equation then reduces to:

(77) $\frac{D\left(L_D\right)}{D_0}=1-\left(1-\frac{1}{?}\right)\frac{\frac{16}{9\sqrt{?}}\frac{L_D}{d}+\left(1-\frac{1}{?}\right)\frac{L_D^2}{L_M^2}}{\frac{16}{9\sqrt{?}}\frac{L_D}{d}+\left(1-\frac{1}{?}\right)\frac{L_D^2}{L_M^2}+\left(1-\frac{1}{?}\right)}.$

(78) FIG. 9 shows the best Pad? fit using the C4(D2O) and C1(D2O) restricted diffusivity at 913 m, and C10 and C1(D2O) at 920 m. C4 and C10 lies in the short diffusion length region, which is sensitive to pore-body size; C1(D2O) lies in the long diffusion length region which is sensitive to tortuosity. Compared with our previous study (Chen et al. 2019) where the largest L.sub.D was 60 ?m, FIG. 9 provides a much more accurate tortuosity since L.sub.D is increased up to 166 ?m.

(79) We note that the gap in L.sub.D between C1 and C10, or between C1 and C4, in FIG. 9 can be closed by using C1 at higher pressure. The higher pressure will decrease the bulk diffusion D.sub.0, which will in term reduce the L.sub.D and close the gap.

(80) According to the Pad? Equation Eq. (5), T.sub.2 and restricted diffusivity are correlated since T.sub.2 is a function of S/V. As such, the Eq. (5) can be converted to a curve on the D-T.sub.2 map. First, the surface relaxivity ? is calculated using d from FIG. 9 and T.sub.2B as such:

(81) $?=(\frac{1}{T_2}-\frac{]}{T_{2B}})\frac{V}{S}=\frac{1}{T_{2S}}\frac{d}{4}$

(82) where T.sub.2 is taken at the mode of region C in FIG. 5, and T.sub.2S is the surface relaxation at the mode of region C. The Pad? fit can then be plotted on the D-T.sub.2 map using the following equation (Zielinski et al., 2010):

(83) 0$\frac{D\left(T_2\right)}{D_0}=1-\left(1-\frac{1}{?}\right)\frac{\frac{4}{9\sqrt{?}}\frac{L_D}{?T_{2S}}+\left(1-\frac{1}{?}\right)\frac{L_D^2}{L_M^2}}{\frac{4}{9\sqrt{?}}\frac{L_D}{?T_{2S}}+\left(1-\frac{1}{?}\right)\frac{L_D^2}{L_M^2}+\left(1-\frac{1}{?}\right)}$

(84) using the best fit values of L.sub.M and T from FIG. 9. In other words, Eq. (14) is the Pad? equation as a function of T.sub.2 (i.e. pore-size) for a fixed L.sub.D, while Eq. (8) is the Pad? equation as a function of L.sub.D for a fixed pore-size d. Both forms of the Pad? equation are equivalent.

(85) FIG. 10 at short diffusion length (L.sub.D=5 ?m) shows the comparison between the Pad? equation and the D-T.sub.2 correlation map for C10. As expected, the Pad? curve shows good agreement with the hydrocarbon-bearing zone in region C, while it starts to deviate from the micritic water signal in region B. This agreement helps validate this technique. Note how using the lower electrical tortuosity (?.sub.E) instead results in a poor fit to the hydrocarbon region, which justifies our invention. In the case of longer diffusion length L.sub.D=13 ?m) in FIG. 10, the electrical tortuosity fails since it yields unphysical ? (i.e. a negative value), which further justifies our invention.

(86) The same analysis can be done for C4(D2O) for the core sample from 913 m, however the bulk T.sub.2B value was not available for C4, which is most likely dominated by dissolved oxygen.

(87) As shown below in the various methods of performing the invention, the Pad? fit to D-T.sub.2 maps is used on NMR logs and core to determine pore-size d, with known L.sub.M and ?.

(88) Uncertainties

(89) As shown above, the Pad? equation can be used to retrieve the pore-body size and tortuosity from NMR restricted diffusivity data. To calculate the error propagation in this regression, the Pad? equation is used as the underlying equation to generate simulated normalized restricted diffusivity. The fixed parameters are chosen in Table 3:

(90) TABLE-US-00003 TABLE 3 Fixed parameter settings in the error analysis. d 5 ?m ? 150 L.sub.M 6 ?m Temperature 30? C. Pressure 1200 psia t.sub.? 15 ? 90 ms

(91) FIG. 11 shows the C10 and C1's restricted diffusivity calculated based on these chosen parameters and its bulk diffusivity. After calculating restricted diffusivity, 1% Gaussian noise is applied to both L.sub.D and the restricted diffusivity. The error in L.sub.D in the experiment mainly comes from the temperature fluctuation during the NMR measurement. The error in the measured normalized restricted diffusivity comes from both the temperature fluctuation and system error in the measurements and inversion algorithm. After getting the noise-added data, the Pad? equation is used to extract the pore-body size and tortuosity information using non-linear curve-fitting (Isqcurvefit) in MATLAB.

(92) Monte Carlo simulations were repeated with 1,000 realizations of the Gaussian noise. Based on 1% Gaussian noise, we calculate an average error ?4.3% for pore-body size and ?1.3% for tortuosity. This indicates that noise has a larger impact on the pore-body-size estimate than on the tortuosity estimate.

(93) The other uncertainty is in selecting the hydrocarbons on the Pad? fit. However as shown in FIG. 12 at 920 m, the best-fit parameters do not change significantly when adding C2, C3, C4, and C5 data at a single measured L.sub.D. More specifically, L.sub.M and d change by ?10%, while ? is more robust and changes by only ?1%.

(94) MICP, BTR, and Permeability

(95) The MICP cumulative saturation C(d) and differential P(d) pore-throat distribution from MICP at 913 m and 920 m are shown in FIG. 13. The total MICP porosity ?.sub.MICP?24 pu agrees with the total NMR porosity ?.sub.T?25 pu for both 913 m and 920 m. The largest d.sub.MICP value (i.e. at the entry pressure), defined at C(d)?0.1%, is d.sub.MICP,entry=0.27 ?m at both 913 m and 920 m, while the mode is d.sub.MICP,mode=0.1 ?m at 913 m and d.sub.MICP,mode=0.07 ?m at 920 m. Note that the micro-pore throat distributions are identical below d.sub.MICP<0.05 ?m for both 913 m and 920 m.

(96) The MICP pore-throat diameters can then be used to determine the body-throat ratio BTR as such:

(97) $\mathrm{BTR}=\frac{d}{d_{\mathrm{MICP},\mathrm{entry}}}\mathrm{or}{\mathrm{BTR}}_{\mathrm{mode}}=\frac{d}{d_{\mathrm{MICP},\mathrm{mode}}}$

(98) where d is the pore-body diameter from the Pad? fit of D/D.sub.0 versus L.sub.D. The results of the BTR are shown in FIG. 14. BTR using d.sub.MICP,entry are lower than BTR.sub.mode using d.sub.MICP,mode. The BTR values are more consistent with expected values for chalks BTR?20 listed in Table 4. The slightly lower BTR value at 913 m is possibly due to the lower d when using C4 (913 m) instead of C10 (920 m), or due to the lower d at 913 m.

(99) TABLE-US-00004 TABLE 4 Body-throat ratios for different rock types according to (Muller-Huber 2016)*, andt he ATR Rock Catalog**. Rock type BTR Microcrystalline dolomite 1-10 * Vuggy dolomite Sucrosic dolomite 10 to 20 * Marly/micritoc limestone Micritic/dolomitic limestone Chalk >20 * Oomoldic limestone Oolitic limestone Sandstones 3-5 **

(100) In cases where MICP data is not available, previously reported values can be used instead, as listed in Table 4 by rock type.

(101) With the region C porosity ?.sub.M, pore-body size d and tortuosity ? all determined from D/D.sub.0 versus L.sub.D, plus the body-throat ratio BTR from d.sub.MICP,entry, the Carman-Kozeny equation predicts the following the permeabilities k:

(102) $k=\frac{{?}_M}{32?}\frac{d^2}{{\mathrm{BTR}}^2}=\frac{0.15}{32?54}?\frac{1.6^2}{6^2}=0.006\mathrm{mD}\left(913m\right)$

(103) $k=\frac{{?}_M}{32?}\frac{d^2}{BT{R}^2}=\frac{0.13}{32?149}?\frac{4.9^2}{18^2}=0.002\mathrm{mD}\left(920m\right)$

(104) These estimates are found to be a factor?? smaller compared to measurements where k.sub.meas=0.035 mD at 913 m, and k.sub.meas=0.017 mD at 920 m. A pre-factor A can be used to match the predictions from measurements, leading to the final expression:

(105) $k=A\frac{{?}_M}{32?}\frac{d^2}{BT{R}^2}\mathrm{where}A?8$

(106) which is used in the methods below to predict permeability on logs and core X.

(107) A small correction in the above formulation which takes the contact angle of mercury on a rough surface into account. On a smooth surface, which is assumed above, the mercury-air contact angle is ?=145?, and the following conversion is used to determine d.sub.MICP from capillary pressure P.sub.c:

(108) $P_c=\frac{4?\cos ?}{d_{MICP}}$

(109) where ? is the mercury-air surface tension. However, for a rough surface has ?=180? (Muller-Huber 2016), resulting in smaller BTR values by a factor 0.81, and therefore A?5 which is closer to unity as predicted by Carman-Kozeny.

(110) Another variation is to use the total porosity ?.sub.T instead of the region C porosity ?.sub.M, which would further reduce the pre-factor to A?2.5.

(111) Validation of Permeability

(112) The validation of the invention using the Carman-Kozeny permeability transform is shown in FIG. 15, for the organic rich chalks from 913 m and 920 m, and for two Austin chalk samples.

(113) For the modified Carman-Kozeny model, a BTR of 5 is assumed for both 913 m and 920 m. The BTR of 5 originates from the Austin chalk sample where MICP pore-throat size (d.sub.throat) using the Swanson method (Swanson, 1981) and pore-body size (d) yield BTR=d/d.sub.throat=5. We do not use d.sub.throat from MICP for 913 and 920 m cores since kerogen/bitumen is highly deformable, which makes the MICP data unreliable, and as shown above elicits an arbitrary pre-factor A. The following estimates are obtained:

(114) $k=\frac{{?}_M}{32?}\frac{d^2}{BT{R}^2}=\frac{0.15}{32?54}?\frac{1.6^2}{5^2}=0.006\mathrm{mD}\left(913m\right)$

(115) $k=\frac{{?}_M}{32?}\frac{d^2}{BT{R}^2}=\frac{0.13}{32?149}?\frac{{4.9}^2}{5^2}=0.026\mathrm{mD}\left(920m\right)$

(116) $k=\frac{{?}_M}{32?}\frac{d^2}{BT{R}^2}=\frac{0.25}{32?7}?\frac{30.6^2}{5^2}=41\mathrm{mD}\left(\mathrm{Austin}\mathrm{Chalk}\right)$

(117) $k=\frac{{?}_M}{32?}\frac{d^2}{{\mathrm{BTR}}^2}=\frac{0.15}{32?9}?\frac{23.3^2}{s^2}=10.8\mathrm{mD}\left(\mathrm{Austin}\mathrm{Chalk}2\right)$

(118) For the Timur-Coates model, the T.sub.2 cutoff for C10 region C is used to separate the BFV and FFV. Two different versions of the SDR model are used; one uses the T.sub.2 logarithmic mean of the SW1 core, and the other uses the T.sub.2 logarithmic mean of C10 region C.

(119) For the sample of Austin chalk, the estimated permeabilities are all near the measured permeabilities. For the low-permeability organic-rich chalk cores from 913 m and 920 m, Timur-Coates and SDR have a large offset compared with measured permeability. The new permeability estimation method uses the modified Carman-Kozeny model with the diffusive tortuosity and pore-body size from NMR restricted diffusion on hydrocarbon-saturated cores with connate water present. This new method estimates the NMR permeability to hydrocarbons at connate water (k.sub.o(S.sub.wc)).

(120) The new method delivers good permeability estimation without empirical parameters, for both the low-permeability organic-rich chalks with their complex pore system and the conventional Austin chalk and Indiana limestone. This new method may also be applied in downhole gradient-based NMR logging for the permeability estimation with the diffusive tortuosity directly measured downhole or from the laboratory.

(121) Core Permeability Using NMR Diffusion-T.sub.2 Measurements

(122) The method requires measuring NMR diffusion-T.sub.2 (D-T.sub.2) and permeability on a selection of cores of different lithologies (i.e. rock types), then using the calibrated parameters (including body-throat-ratio BTR and cementation exponent m) to determine permeability of core X in a large case study, without having to measure a (more time consuming) light hydrocarbon (e.g. C1) for each core X.

(123) Method: 1. Obtain cores from a selection of depths to cover all lithologies (i.e. rock types). 2. Measure D-T.sub.2 on cores to determine movable fluid porosity (?.sub.M), i.e. porosity above T.sub.2>T.sub.2cutoff, and restricted diffusion (D/D.sub.0), where D.sub.0 is bulk diffusion versus diffusion length

(124) 0$(L_D=\sqrt{D_0{t}_{?}},$ where t.sub.? is the diffusion encoding time) using a light hydrocarbon (e.g. C1) and a heavier hydrocarbon (e.g. C4 or C10) to cover a range of L.sub.D. 3. Fit D/D.sub.0 versus L.sub.D using Pad? fit to determine pore body size (d), heterogeneity length scale (L.sub.M), and tortuosity (?) for each core. 4. Determine the cementation exponent (m) for each lithology using ?=?.sub.M.sup.1-m. 5. Determine heterogeneity length scale (L.sub.M) for each lithology. 6. Measure core permeability (k) and determine the body-throat ratio BTR for each lithology from the Carman-Kozeny equation

(125) $k=\frac{{?}_M{d}^2}{32?BT{R}^2}$ 7. Measure D-T.sub.2 on core X to determine movable fluid porosity (?.sub.M,X), i.e. porosity above T.sub.2>T.sub.2cutoff, and D.sub.X/D.sub.0 versus T.sub.2 using a heavier hydrocarbon (e.g. C4 or C10) to cover a short L.sub.D. 8. Determine tortuosity log (?.sub.X) using

(126) ${?}_X={?}_{M,X}^{1-m},$ where m is known for each lithology from step 4. 9. Fit D.sub.X/D.sub.0 versus T.sub.2 using Pad? fit to determine pore-body size (d.sub.X), with ?.sub.X known from step 8, and L.sub.M fixed for each lithology from step 5. 10. Determine permeability (k.sub.X) for core X from

(127) $k_X=\frac{{?}_{M,X}{d}_X^2}{32{?}_X{\mathrm{BTR}}^2},$ with BTR known for each lithology from step 6.

(128) Variations: 1. Use total porosity (?.sub.T) instead of movable fluid porosity (?.sub.M). 2. Use an amplitude factor (A) in the Carman-Kozeny relation

(129) $k=A\frac{{?}_M{d}^2}{32?BT{R}^2}$ 3. Fix tortuosity (?) instead of cementation exponent (m) for each lithology. 4. Do not use heterogeneity length-scale (L.sub.M) in Pad? fit, i.e. assume L.sub.D/L.sub.M=0. 5. Use another fitting function for D/D.sub.0 versus L.sub.D besides the Pad? fit. 6. Use another permeability prediction besides Carman-Kozeny relation. 7. Fix BTR based on known values for each lithology. 8. Determine BTR=d/d.sub.throat from MICP d.sub.throat for each lithology. 9. Use different fluids with a large D.sub.0 (i.e., not C1) and small D.sub.0 (i.e., not C4 and not C10). 10. Make no distinction in lithology. 11. Deuterate the connate water with heavy-water or dope it with paramagnetics to remove the water signal. 12. Measure the core at 100% hydrocarbon saturation without connate water.

(130) Downhole Permeability Log Using NMR Diffusion-T.sub.2 Measurements

(131) The method requires measuring NMR diffusion-T.sub.2 (D-T.sub.2) and permeability on a selection of cores of different rock types, then using the core calibrated parameters (including body-throat-ratio BTR and cementation exponent m) to interpret D-T.sub.2 logs to obtain a permeability log.

(132) More specifically, the permeability is measured in the laboratory and used with the permeability equation to determine a value of BTR for that lithology. That BTR and m for each lithology are then used in the downhole wireline NMR diffusion logging with a gradient-based NMR logging tool to generate a continuous downhole permeability log. Examples of gradient-based downhole NMR logging tool include Baker-Hughes MREX, Schlumberger's MR Scanner and Halliburton's MRIL.

(133) Method: 1. Obtain cores from a selection of depths to cover all lithologies (i.e. rock types). 2. Measure D-T.sub.2 on cores to determine movable fluid porosity (?.sub.M), i.e. porosity above T.sub.2>T.sub.2cutoff, and restricted diffusion (D/D.sub.0), where D.sub.0 is bulk diffusion versus diffusion length

(134) $(L_D=\sqrt{D_0{t}_{?}},$ where t.sub.? is the diffusion encoding time) using a light hydrocarbon (e.g. C1) and a heavier hydrocarbon (e.g. C4 or C10) to cover a range of L.sub.D. 3. Fit D/D.sub.0 versus L.sub.D using Pad? fit to determine pore body size (d), heterogeneity length scale (L.sub.M), and tortuosity (?) for each core. 4. Determine the cementation exponent (m) for each lithology using ?=?.sub.M.sup.1-m. 5. Determine heterogeneity length scale (L.sub.M) for each lithology. 6. Measure core permeability (k) and determine the body-throat ratio BTR for each lithology from the Carman-Kozeny equation

(135) $k=\frac{{?}_M{d}^2}{32?BT{R}^2}$ 7. Measure D-T.sub.2 log to determine movable fluid porosity (?.sub.M, log), i.e. porosity above T.sub.2>T.sub.2cutoff, and D.sub.log/D.sub.0 versus T.sub.2. 8. Determine tortuosity log (?.sub.log) using

(136) ${?}_{log}={?}_{M,\log }^{1-m},$ where m is known for each lithology from step 4. 9. Fit D.sub.log/D.sub.0 versus T.sub.2 using Pad? fit to determine pore-body size (d.sub.log), with ?.sub.log known from step 8, and L.sub.M fixed for each lithology from step 5. 10. Determine downhole permeability (k.sub.log) from

(137) $k_{\log }=\frac{{?}_{M,\log }{d}_{\log }^2}{32{?}_{\log }{\mathrm{BTR}}^2},$ with BTR known for each lithology from step 6.

(138) Variations: 1. Use total porosity (?.sub.T) instead of movable fluid porosity (?.sub.M). 2. Use an amplitude factor (A) in the Carman-Kozeny relation

(139) $k=A\frac{{?}_M{d}^2}{32?BT{R}^2}$ 3. Fix tortuosity (?) instead of cementation exponent (m) for each lithology. 4. Do not use heterogeneity length-scale (L.sub.M) in Pad? fit, i.e. assume L.sub.D/L.sub.M=0. 5. Use another fitting function for D/D.sub.0 versus L.sub.D besides the Pad? fit. 6. Use another permeability prediction besides Carman-Kozeny relation. 7. Fix BTR based on known values for each lithology. 8. Determine BTR=d/d.sub.throat from MICP d.sub.throat for each lithology. 9. Use different fluids with a large D.sub.0 (i.e., not C1) and small D.sub.0 (i.e., not C4 and not C10). 10. Make no distinction in lithology. 11. Deuterate the connate water with heavy-water or dope it with paramagnetics to remove the water signal. 12. Measure the core at 100% hydrocarbon saturation without connate water.

(140) New Permeability Equation Extending BTR

(141) New permeability equations may be based on empirical findings of the data, or by modifying the Carman-Kozeny (or other) equation based on physical models. One such modification based on physical models incorporates work by Monicard 1980, see FIG. 16, where instead of parallel cylinders, the flow paths are comprised of a sequence of cone frustrums of pore-body size d and pore-throat size d.sub.throat=d/BTR. The BTR value characterizes the ratio of the widest to narrowest sections of the cones.

(142) The modified Carman-Kozeny equation is then:

(143) 0$k=A\frac{{?}_M}{32?}\frac{d^2}{BT{R}^2}M\mathrm{where}M=\frac{9BT{R}^3}{{\left(BT{R}^2+BTR+1\right)}^2}$

(144) The extra term M is equal to unity M=1 at BTR=1, tends towards M?9 BTR for BTR<<1, and tends towards M?9/BTR for BTR>>1, as shown in FIG. 17.

(145) This equation provides a new method for determining BTR in the laboratory. The permeability k of the core is measured by routine core analysis (either steady-state or unsteady-state methods), ?.sub.M, ?, and d are measured by NMR diffusion, set A=1, and this equation is then solved for the BTR of the mobile hydrocarbons. This new method of determining BTR is non-destructive (no Hg injection) and more accurately determines the BTR of the portion of the pore size distribution where the mobile hydrocarbons reside. The BTR determined using this method can be applied to the downhole wireline NMR log for determining a continuous permeability log of the mobile hydrocarbons.

(146) Variations on the Methodology

(147) Saturation State of the Core

(148) The first option for core saturation is to measure the core at irreducible water saturation S.sub.wirr, which is the case shown in FIG. 5 and all subsequent figures. In this case, the connate water is present in Region B, and the tortuosity of the hydrocarbon phase is measured in Region C. In conventional reservoirs, having the core at S.sub.wirr resembles reservoir conditions at heights well above the free water level. In other words, having the core at S.sub.wirr resembles downhole conditions in the oil-bearing zone, and therefore allows for better comparison with permeability from log data.

(149) Another reason for having the core at S.sub.wirr is that in the case of chalks, the micritic porosity has a large tortuosity which does not contribute to permeability. As such, the tortuosity from the micritic porosity should not be counted when deriving a permeability transform which uses tortuosity as an input. The best way to avoid this is to have the core at S.sub.wirr and to measure tortuosity of the movable hydrocarbon (which is immiscible with water), as illustrated in FIG. 4. In such cases, the tortuosity of the hydrocarbon phase in Region C will correspond to the tortuosity used in a permeability transform. Note that measuring tortuosity of the hydrocarbon phase at S.sub.wirr is possible with NMR, but it is not possible with electrical resistivity.

(150) The second option for core saturation is to measure the core 100% saturated with a single fluid, e.g. 100% brine saturated or 100% hydrocarbon saturated. This option is appropriate when there is no micritic porosity, as in the case for sandstones, limestones, and dolomites. In such cases, the tortuosity of the 100% saturating phase corresponds to the tortuosity in a permeability transform, therefore the core does not need to be at S.sub.wirr.

(151) FIG. 18 below shows the results of 100% H2O saturated data, and 100% C1 saturated data after drying, on an Indiana limestone core.

(152) Using the above Pad? fit, the following calculation is made:

(153) $k=\frac{{?}_M}{32?}\frac{d^2}{BT{R}^2}=\frac{0.17}{32?3}?\frac{21.6^2}{5^2}=32\mathrm{mD}\left(\mathrm{Indiana}\mathrm{limestone}\right)$

(154) which shows excellent agreement with the range of measured values for this block of Indiana limestone ranging between 30-50 mD.

(155) Determining the Diffusion Coefficient

(156) The diffusion coefficient D used to plot restriction D/D.sub.0 versus L.sub.D in FIG. 12 was determined by the 2D peak (i.e. mode) of the D-T.sub.2 map, hereby called D.sub.2Dpeak. Another way to determine D is to fit the initial slope of the magnetization decay due to diffusion encoding, as shown in the example in FIG. 19 for an Austin Chalk sample at S.sub.wirr with C10 and D.sub.2O connate brine. In this example, the diffusion coefficient from the initial slope D.sub.slope=0.781 ?m.sup.2/ms is ?40% lower than D.sub.2Dpeak=1.237 ?m.sup.2/ms at the 2D peak of the D-T.sub.2 map.

(157) Yet another way to determine the diffusion coefficient is at the peak of the 1D diffusion distribution D.sub.1Dpeak=1.102 ?m.sup.2/ms, or the log-mean value of the 1D diffusion D.sub.Lm=0.727 ?m.sup.2/ms. The log-mean method is comparable to the slope method, while the 1D and 2D peak methods yield ?30-40% larger diffusion coefficients.

(158) One big caveat with the slope, 1D peak and log-mean methods for cores at S.sub.wirr is that the cores must have deuterated (D.sub.2O) connate brine, otherwise these techniques will yield an average diffusion coefficient from the H.sub.2O connate brine and the hydrocarbon phase, which is incorrect. The 2D peak technique does not suffer from this issue since a T.sub.2 cutoff is applied to separate hydrocarbons from water before determining D.sub.2Dpeak.

(159) An example of the effects of D.sub.slope versus D.sub.2Dpeak to the Pad? fit is shown in FIG. 20. The Pad? fit using D.sub.2Dpeak indicates a reasonable tortuosity of ?=7 but a large pore-size of d=30.6 ?m. On the other hand, using D.sub.slope indicates a reasonable pore-size of d=7.5 ?m but a large tortuosity of ?=20. Further investigations are required to determine which measure of D is more accurate, however note that the Carman-Kozeny permeability transform is more susceptible to inaccuracies in pore-size (k?d.sup.2) than to tortuosity (k?1/?).

(160) Another caveat with the slope technique is that D.sub.slope represents the average diffusion coefficient over all pore sizes, without the option of selecting pore types based on T.sub.2. This presents an issue in cores which contain a significant amount of vugs, an example of which is shown in FIG. 21 for an Indiana limestone core saturated with brine. The large peak at long T.sub.2>750 ms originates from vugs with almost no surface relaxation and no restricted diffusion. Furthermore, vugs do not contribute to the NMR permeability transform, therefore only signal below T.sub.2<750 ms should be considered (Chang et al, 1994, Hidajat et al. 2004).

(161) Using the D.sub.slope in such cases results in an artificially low tortuosity since the vugs dominate the signal. This accounts for the lower tortuosity from D.sub.slope compared to electrical tortuosity in Indiana Limestone, as was previously reported (Yang et al., 2019). The solution is to use D.sub.2Dpeak in the region T.sub.2<750 ms instead. Note that D.sub.1Dpeak and D.sub.LM would also fail in such cases.

(162) Determining the Tortuosity

(163) Another method for determining tortuosity is not to use the Pad? model (or any other model), but to simply determine the tortuosity from the average of the C1(D2O) data at large L.sub.D. An example of this technique is shown in FIG., where it is clear that the average tortuosity is smaller than the Pad? fit results from FIG. 20. In the case of D.sub.2Dpeak the average tortuosity is ?30% smaller compared to the Pad? tortuosity, while in the case of D.sub.slope the average tortuosity is ?60% smaller compared to the Pad? tortuosity. Note that the y-axis in FIG. 22 is on a linear scale instead of a log scale.

(164) The advantage of the average tortuosity is that no model (such as the Pad? model) is required to fit the data. The disadvantage of the average tortuosity method is that the pore-size cannot be determined since there is no information used at short L.sub.D. A potential work around is to use a separate fit to the initial slope of D/D.sub.0 versus L.sub.D to get the average pore-size, however as shown in FIG. 22, the decay at the shortest L.sub.D is already too large to fit a linear decay, thereby yielding an inaccurate pore-size.

(165) Improving Accuracy at Short Diffusion Lengths

(166) Another variation of the above methodology is to obtain more accuracy at short L.sub.D by reducing D.sub.0 of the bulk alkane, thereby improving the accuracy of the pore-size estimation. This can be accomplished by increasing the carbon number above C10 and/or reducing the temperature, thereby reducing D.sub.0 and L.sub.D.

(167) TABLE-US-00005 TABLE 5 List of viscosities (Rossini et al. 1953), bulk diffusion coefficients (Tofts et al. 2000), and diffusion lengths (at the shortest encoding time) of various n-alkanes, and water. T P ? D.sub.0 t.sub.? L.sub.D Molecule (C) (bar) (cP) (?m.sup.2/ms) (ms) (?m) H2O 23 1 1.00 2.15 15 5.68 H2O 4 1 1.52 1.28 15 4.38 C4 23 2.2 0.18 7.54 15 10.64 C5 20 1 0.24 4.34 15 8.07 C6 23 1 0.31 3.70 15 7.45 C7 23 1 0.42 2.99 15 6.70 C8 20 1 0.55 2.14 15 5.67 C9 20 1 0.72 1.63 15 4.94 C10 20 1 0.93 1.26 15 4.35 C10 (*) -10 1 1.39 0.76 15 3.37 C11 20 1 1.19 1.01 15 3.89 C12 20 1 1.51 0.78 15 3.42 C13 20 1 1.89 0.63 15 3.07 C14 20 1 2.34 0.49 15 2.71 C15 20 1 2.87 0.40 15 2.45 C16 20 1 3.48 0.34 15 2.26 (*) viscosity extrapolated from Lee et al. 1965, assuming D.sub.0 ? T/?.

(168) As listed in Table 5, increasing the carbon number from C10 to C16 reduces D.sub.0 by a factor?4, and therefore reduces L.sub.D by a factor?2 from L.sub.D?4 ?m to L.sub.D?2 ?m. This will significantly improve the pore-size estimation for pares less than d<5 ?m.

(169) Another way to reduce D.sub.0 and L.sub.D is to lower the temperature of the measurement. As listed in Table 5, if water is used as the saturating fluid, D.sub.0 decreases by ?50% and L.sub.D decreases by ?25% when going from ambient temperature to 4? C. If decane is used as the saturating fluid, D.sub.0 decreases by ?40% and L.sub.D decreases by ?20% when going from ambient temperature to ?10? C.

(170) An extension of this technique is to reduce the temperature below T<0? C. such that the water signal disappears, i.e. T.sub.2<t.sub.E for ice, and one is left with only the hydrocarbon NMR signal. This will replace the need for D.sub.2O exchange with water in the above methodology, and it will reduce L.sub.D of the hydrocarbon for better determination of pore-size. One caveat is that the water may freeze below T<0? C. in nano-pores due to the Gibbs-Thomson effect, in which case lower temperatures are required; this phenomenon can (in principle) also be used to determine the nano-pore size distribution in porous media, a.k.a. NMR cryoporometry (Weber 2010).

(171) Fluids with Large Diffusivity

(172) A fluid with large diffusivity D.sub.0 (i.e. large diffusion length L.sub.D) is required for determining the tortuosity in the Pad? fit. For the fluid candidate which has large bulk diffusivity, CH.sub.4 is a good candidate, when using .sup.1H-NMR. Even though H.sub.2 and HD have larger bulk diffusivities than methane, they have lower HI. FIG. 23 shows the hydrogen index (HI) comparison between water, methane, H.sub.2, and HD. We obtain the molar density of the fluids from the REFPROP database. The temperature for the fluid properties is fixed at 30? C. The pressure ranges from 1 psia to 2,000 psia. HD's HI is too low. H.sub.2's HI is relatively higher than HD. And for bulk H.sub.2, T.sub.2?15 ms at 147 psi and room temperature. For the high-pressures, the T.sub.2 should be higher than 15 ms, because T.sub.2 is linearly dependent on the density (Abragam, 1961).

(173) Workflow to Optimize Pulse Parameters

(174) The workflow shown below consists of a quick measurement with default pulse parameters, followed by a full measurement with optimized parameters: 1. Measure a quick D-T.sub.2 using the default pulse parameters, which comprises a minimum number of scans (NS), a small number of gradient steps (NG), and the minimum diffusion time (t.sub.?). The gradient list in the range G=0.fwdarw.G.sub.max is set based on NG linearly-spaced steps, where G.sub.max is the maximum gradient strength. 2. Optimize G.sub.max based on diffusion decay M(b) by either (1) increasing G.sub.max if more decay is required in M(b) due to slowly decaying components, (2) decreasing G.sub.max if there is insufficient data at short b due to fast decaying components. An algorithm can be used to forward model M(b) based on D-T.sub.2 from the default parameters, from which G.sub.max can be optimized. 3. Measure a full D-T.sub.2 using a larger number of scans (NS) for higher SNR, and a gradient list in the range G=0.fwdarw.G.sub.max using the optimized G.sub.max, with a larger number of linearly-spaced gradient steps (NG) to obtain more data at short b. 4. Measure full D-T.sub.2 at larger diffusion-times (t.sub.?) in order to probe larger diffusion-lengths (L.sub.D). For each t.sub.?, G.sub.max is now optimized by keeping

(175) $b_{\max }=q_{\max }^2\left(t_{?}-t_{?}/3+?\right)$ roughly constant, where

(176) $q_{\max }=?t_{?}{G}_{\max }.$ In other words, the product t.sub.?G.sub.max.sup.2 is kept roughly constant for each t.sub.? measurement.

(177) An example of the diffusion decay M(b) using the default parameters is shown in FIG. 24, as well as the default parameters used. The parameters are defined in FIG. 1. NG is defined as the number of gradient steps, and NS is the number of scans which according to Table 1 has to be a multiple of 16 due to phase cycling. The gradient list in the range G=0.fwdarw.G.sub.max is set based on NG linearly-spaced steps, where G.sub.max is the maximum gradient strength. The minimum NS=16 is used in order to save time. The minimum t.sub.? is used in order to capture the most signal (i.e. the shortest dead time for encoding diffusion).

(178) The above core example has both irreducible water represented by the long tail in the M(b) decay with small D, and methane signal represented by the fast-decaying portion in M(b) with large D. Note that the linear spacing in the gradient list G means that there are more points clustered at small b compared to large b since b?G.sup.2. This is intentional since typically more resolution is required at small b.

(179) The quick test in FIG. 24 shows insufficient decay in M(b); from ?2.7 pu to ?0.7 pu (i.e. only a factor?4 decay), which will limit the resolution in the slow decaying component. As such, G.sub.max is increased from 12 G/cm to 20 G/cm in order to increase the measured decay in M(b). The optimized data are shown in FIG. 25, which now decays from ?2.7 pu to ?0.4 pu, i.e. by a factor?7 decay. Increasing G.sub.max further (i.e. to the instrument maximum of 45 G/cm) is not optimal since a sufficient number of points in the gradient list are required at short b to capture the initial fast-decaying component in M(b). For instance, the initial decay in FIG. 25 is captured by the first ?16 data points alone. An algorithm can be used to forward model M(b) based on D-T.sub.2 from the default parameters, from which G.sub.max can be optimized.

(180) The difference between the default and optimal parameters on the diffusion projection is shown in FIG. 26. The insufficient decay in M(b) using the default parameters results in a long tail in P(D) at small D (i.e. the irreducible water), which is remedied when using the optimal parameters. Note also how the peak in P(D) at large D (i.e. the methane) is slightly sharper for the optimized parameters, which is a result of having higher SNR and more data points in the initial decay of M(b), that is ?16 data points versus ?9.

(181) The difference between the default and optimal parameters on the on D-T.sub.2 maps is shown in FIG. 27, from which the same conclusions can be reached as from the diffusion projection in FIG. 26.

(182) Downhole NMR Logging Applications

(183) The modified Carman-Kozeny permeability estimation method using hydrocarbon diffusive tortuosity from D-T.sub.2 measurements can be applied to downhole NMR logging applications. There are several practical considerations listed below: 1. This invention needs fluids with different bulk diffusivity to cover a wide range of the diffusion length (L.sub.D). In the downhole NMR logging application, if there are gas and liquid, we can distinguish them using a D-T.sub.2 map and estimate their restricted diffusivity (Minh et al. 2015, Vinegar et al. 2020). 2. This invention has been demonstrated on (organic-rich) chalk and other carbonate formations such as limestones and dolomites. For maximum accuracy, the BTRs should be determined in the laboratory for different formations using this invention and those BTRs should be used in the downhole NMR logging.

(184) This invention may require relatively high signal-to-noise ratio (SNR) in D-T.sub.2 measurements. Indeed, in some cases the downhole NMR gradient-based logging tools may not get such high-quality data as obtained in the laboratory (Minh et al. 2015, Vinegar et al. 2020). In the NMR downhole measurements, taking an average of adjacent depths may be necessary to increase SNR. Holding the logging tool stationary to obtain a large number of averages is another way to increase SNR.

(185) The present invention has been described using detailed descriptions of embodiments thereof that are provided by way of example and are not intended to limit the scope of the invention. The described embodiments comprise different features, not all of which are required in all embodiments of the invention. Some embodiments of the present invention utilize only some of the features or possible combinations of the features. Variations of embodiments of the present invention that are described and embodiments of the present invention comprising different combinations of features noted in the described embodiments will occur to persons skilled in the art.

REFERENCES

(186) ATR Rock Catalog Archie, G. E., 2003: The Electrical Resistivity Log as an Aid in Determining Some Reservoir Characteristics. SPE Reprint Series, 146, 9-16. Bernab?, Y.; Mok, U.; Evans, B., 2006: A Note on the Oscillating Flow Permeability. International Journal of Rock Mechanics and Mining Sciences, 43, 311-316. Carey, G. R.; McBean, E. A.; Feenstra, S., 2016: Estimating Tortuosity Coefficients Based on Hydraulic Conductivity. Groundwater, 54, 476-487. Carman, P. C., 1997: Fluid Flow Through Granular Beds, Process Safety and Environmental Protection: Transactions of the Institution of Chemical Engineers, Part B., 75, S32-S48. Chang, D., Vinegar, H. J., Morriss, C., Straley, C., 1994: Effective Porosity, Producible Fluid and Permeability in Carbonates from NMR Logging, Society of Petrophysicists and Well-Log Analysts, SPWLA-1994-A Chen, Z.; Singer, P. M.; Kuang, J.; Vargas, F. M.; Hirasaki, G. J., 2017: Effects of Bitumen Extraction on the 2D NMR Response of Saturated Kerogen Isolates. Petrophysics, 58, 470-484. Chen, Z., Singer, P. M., Wang, X., Vinegar, H. J., Nguyen, S. V., Hirasaki, G. J., 2019: NMR Evaluation of Light-Hydrocarbon Composition, Pore Size, and Tortuosity in Organic-Rich Chalks, Petrophysics 60 (06), 771-797 Chen, Z., Singer, P. M., Wang, X., Hirasaki, G. J., Vinegar, H. J., 2019: Evaluation of Light Hydrocarbon Composition, Pore Size, and Tortuosity in Organic-Rich Chalks using NMR Core Analysis and Logging, Society of Petrophysicists and Well-Log Analysts, SPWLA-2019-K Chen, Z., Wang, X., Jian, G., Zhang, L., Dong, P., Singer, P. M., & Hirasaki, G. J., 2020: Fast Permeability Estimation of an Unconventional Formation Core by Transient-Pressure History Matching, Society Petroleum Engineers Journal Coates G. R., Miller, M., Gillen, M., Henderson, C., 1991: The MRIL in Conoco 33-1: An Investigation of a New Magnetic Resonance Imaging Log, Society of Petrophysicists and Well-Log Analysts, SPWLA-1991-DD Cotts, R. M., Hoch, M. J. R., Sun, T., Markert, J. T.; 1989: Pulsed Field Gradient Stimulated Echo Methods for Improved NMR Diffusion Measurements in Heterogeneous Systems. J. Magn. Reason. 83(2), 252-266 Dullien, F. A. L., 1979: Porous Media Fluid Transport and Pore Structure. Academic Press Epstein, N., 1989: On Tortuosity and the Tortuosity Factor in Flow and Diffusion Through Porous Media, Chemical Engineering Science., 44, 777-779. Glover, P. W. J., 2016: Archie's LawA Reappraisal. Solid Earth., 7, 1157-1169. Hidajat, I., Mohanty, K. K., Flaum, M., Hirasaki, G. J., 2004, Study of Vuggy Carbonates Using NMR and X-Ray CT Scanning, SPE Reservoir Evaluation & Engineering, SPE 88995 Hurlimann, M. D.; Helmer, K. G.; Latour, L. L.; Sotak, C. H., 1994: Restricted Diffusion in Sedimentary Rocks. Determination of Surface-Area-to-Volume Ratio and Surface Relaxivity, Journal of Magnetic Resonance, Series A., 111, 169-178. Kenyon, W. E., Day, P. I., Straley, C., Willemsen, J. F., 1988: A Three-part Study of NMR Longitudinal Relaxation Properties of Water-saturated Sandstones, Society Petroleum Engineers, SPE-15643-PA Latour, L. L.; Mitra, P. P.; Kleinberg, R. L.; Sotak, C. H., 1993: Time-Dependent Diffusion Coefficient of Fluids in Porous Media as a Probe of Surface-to-Volume Ratio. Journal of Magnetic ResonanceSeries A., 101, 342-346. Lee, A. L., Ellington, R. T., 1965, Viscosity of n-Decane in the Liquid Phase, Journal of Chemical and Engineering Data, 346-348 Lo, S. W.; Hirasaki, G. J.; House, W. V.; Kobayashi, R., 2002: Mixing Rules and Correlations of NMR Relaxation Time with Viscosity, Diffusivity, and Gas/Oil Ratio of Methane/Hydrocarbon Mixtures. SPE Journal., 7, 24-34. Minh, C. C., Crary, S., Singer, P. M., Valori, A., Bachman, N., Hursan, G. G., Ma, S., Belowi A., Kraishan G., Kraishan, G. (2015). Determination of Wettability from Magnetic Resonance Relaxation and Diffusion Measurements on Fresh-State Cores. SPWLA 56th Annual Logging Symposium. Mitchell, J.; Gladden, L. F.; Chandrasekera, T. C.; Fordham, E. J., 2014: Low-Field Permanent Magnets for Industrial Process and Quality Control. Progress in Nuclear Magnetic Resonance Spectroscopy., 76, 1-60. Mitra, P. P.; Sen, P. N., 1992: Effects of Microgeometry and Surface Relaxation on NMR Pulsed-Field-Gradient Experiments: Simple Pore Geometries. Physical Review B., 45, 143-156. Monicard, R. P., 1980: Properties of Reservoir Rocks: Core Analysis. Springer, Netherlands. Muller-Huber, E., Schon, J., Borner, F., 2016 A Pore Body-Pore Throat-Based Capillary Approach for NMR Interpretation in Carbonate Rocks using the Coates Equation, Society of Petrophysicists and Well-Log Analysts Rossinni, F. D., Pitzer, K. S., Arnett, R. L., Braun, R. M., Pimente, G. C., 1953 Selected Values of Physical and Thermodynamic Properties of Hydrocarbons and Related Compounds, Tech. Rep. Project 44, American Petroleum Institute, Pittsburgh, PA Sander, R.; Pan, Z.; Connell, L. D., 2017: Laboratory Measurement of Low Permeability Unconventional Gas Reservoir Rocks: A Review of Experimental Methods. Journal of Natural Gas Science and Engineering., 37, 248-279. Singer, P.; Chen, Z.; Hirasaki, G., 2016: Fluid Typing and Pore Size in Organic Shale Using 2D NMR in Saturated Kerogen Isolates. Petrophysics., 57, 604-619. Singer, J. M., Johnston, L., Kleinberg, R. L., Flaum, C., 1997, Fast NMR Logging for Bound Fluid and Permeability, Society of Petrophysicists and Well-Log Analysts, SPWLA-1997-YY Souza, A., Carneiro, G., Zielinski, L., Polinski, R., Schwartz, L., Hurlimann, M. D., Boyd, A., Rios, E. H., Camilo dos Santos, B. C., Trevizan, W. A., Machado, F. V., Azeredo, R. B. V., 2013: Permeability Prediction Improvement using 2D NMR Diffusion-T.sub.2 Maps, Society of Petrophysicists and Well-Log Analysts, SPWLA-2013-U Straley, C., Rossini, D., Vinegar, H. J., Tutunjian, P., Morriss, C., 1997: Core Analysis by Low-Field NMR, Society of Petrophysicists and Well-Log Analysts, SPWLA-1997-v38n2a3 Swanson, B. F., 1981: A Simple Correlation Between Permeabilities and Mercury Capillary Pressures. Society of Petroleum Engineers. Journal of Petroleum Technology., 33, 2498-2504. Tanner, J. E., 1970: Use of the Stimulated Echo in NMR Diffusion Studies. The Journal of Chemical Physics., 52, 2523-2526. Timur, R., 1969: Pulsed Nuclear Magnetic Resonance Studies of Porosity, Movable Fluid, and Permeability of Sandstones, Society Petroleum Engineers, SPE-2045-PA Tofts, P. S., Lloyd, D., Clark, C. A., Barker, G. J., Parker, G. J. M., McConville, P., Baldock, C., Pope, J. M., 2000 Test Liquids for Quantitative MRI Measurements of Self-Diffusion Coefficient In Vivo, Magn. Reson. Med. 43, 368-374 Venkataramanan, L.; Song, Y. Q.; H?rlimann, M. D., 2002: Solving Fredholm integrals of the first kind with tensor product structure in 2 and 2.5 dimensions. IEEE Transactions on Signal Processing., 50, 1017-1026. Vinegar, E. G.; Rosenberg, Y. O.; Reznick, I.; Gordin, Y.; Singer, Philip M. S.; Wang, X.; Chen, Z.; Nguyen, S. V.; Li, W.; Bradley, T.; Hirasaki, G. J.; Lake, L. W.; Feinstein, S.; Hatzor, Y. H.; Vinegar, H. J., 2020. What Happens to the Petrophysical Properties of a Dual-Porosity Organic-Rich Chalk During Early-Stage Organic Maturation? Society of Petrophysicists and Well-Log Analysts. Wang, X., Singer, P. M., Chen, Z., Hirasaki, G. J., Vinegar, H. J., 2020: A New Method of Estimating Tortuosity and Pore Size in Unconventional Formations using NMR Restricted Diffusion Measurements, Society of Petrophysicists and Well-Log Analysts Webber, J. B. W., 2010 Studies of Nano-Structured Liquids in Confined Geometries and at Surfaces, Progress in Nuclear Magnetic Resonance Spectroscopy 56 78-93 Wyllie, M. R. J.; Rose, W. D., 1950: Some Theoretical Considerations Related To The Quantitative Evaluation Of The Physical Characteristics Of Reservoir Rock From Electrical Log Data, Journal of Petroleum Technology., 2, 105-118 Yang, K.; Li, M.; Ling, N. N. A.; May, E. F.; Connolly, P. R. J.; Esteban, L.; Clennell, M. B.; Mahmoud, M.; El-Husseiny, A.; Adebayo, A. R.; Elsayed, M. M.; Johns, M. L., 2019: Quantitative Tortuosity Measurements of Carbonate Rocks Using Pulsed Field Gradient NMR, Transport in Porous Media, 130, 847-865. Zecca, M.; Vogt, S. J.; Connolly, P. R. J.; May, E. F.; Johns, M. L., 2018: NMR Measurements of Tortuosity in Partially Saturated Porous Media. Transport in Porous Media, 125, 271-288. Zielinski, L.; Ramamoorthy, R.; Minh, C. C.; Al Daghar, K. A.; Sayed, R. H.; Abdelaal, A. F., 2010: Restricted Diffusion Effects in Saturation Estimates From 2D Diffusion-Relaxation NMR Maps. Society of Petroleum Engineers.

ABBREVIATIONS

(187) ASREC: as-received core C10: decane-saturated core C4: butane-saturated core C1(D2O): methane-saturated core with D2O D2O: deuterated core NMR: Nuclear Magnetic Resonance PFG: pulse-field gradient RF: radio frequency CPMG: Carr-Purcell-Meiboom-Gill sequence PTFE: polytetrafluoroethylene PEEK: polyether ether ketone HI: hydrogen index pu1: porosity unit from 1D T.sub.2 assuming HI=1 pu1.sub.D: porosity unit from D-T.sub.2 assuming HI=1 SW1: 100% brine saturated

SYMBOLS

(188) ?.sub.e: electrical tortuosity F.sub.R: resistivity formation factor ?: porosity m: Archie cementation exponent L.sub.e: true travel length for molecule L: direct geometry ?: diffusive tortuosity D.sub.0: free self-diffusion of bulk fluid D.sub.?: restricted self-diffusion in the tortuosity limit D: restricted diffusivity D.sub.slope: diffusion coefficient from initial slope of diffusion decay D.sub.1Dpeak: diffusion coefficient from peak of D distribution D.sub.2Dpeak: diffusion coefficient from peak of D-T.sub.2 map D.sub.LM: diffusion coefficient from log-mean of D distribution ?.sub.?: diffusion time S/V: surface-to-volume ratio L.sub.D: diffusion length L.sub.M: heterogeneity length-scale d: cylindrical pore-body diameter d.sub.MICP: cylindrical pore-body diameter T.sub.1: longitudinal relaxation time T.sub.1LM: log mean longitudinal relaxation time T.sub.2: transverse relaxation time T.sub.2B: transverse relaxation time of bulk fluid T.sub.2s: surface transverse relaxation time t.sub.?: duration of the gradient encoding pulses G: gradient in trapezoidal gradient encoding pulses ?: ramp-time t.sub.D: dead time ?.sub.1,2,3,a: RF phases q: wavevector for diffusion encoding ?/2?: gyro-magnetic ratio P: probability density function ?.sub.DT.sub.2: total porosity from D-T.sub.2 map ?: logarithmic bin spacings of the 2D map B: boosting factor t.sub.E: echo space ?: surface relaxivity