CN116507913A - Direct flow simulation of calibrated rock samples - Google Patents

Abstract

Disclosed herein is a method for calibrating a direct flow simulation of a rock sample, the method involving providing a 3D image of the rock sample, and generating a segmented structural image of the rock sample from the 3D image by selecting voxels to represent void space or solid material. Fluid flow is simulated on the segmented structure image using direct flow simulation. A 3D spatially resolved fluid velocity map of one or more fluid phases is generated at a pore scale resolution using pulsed field gradient magnetic resonance imaging. The simulated fluid flow and the 3D spatially resolved fluid velocity map are compared to calibrate a direct flow simulation on the rock sample.

Description

Direct flow simulation of calibrated rock samples

Technical Field

The present disclosure relates to the field of direct flow simulation of rock samples. In particular, the method calibrates direct flow simulation with 3D images generated by digital petrophysical enhanced by flow magnetic resonance imaging (flow MRI).

Background

Analysis and imaging techniques that can help visualize fluid flow in the pore space of sedimentary rock are of great interest to the oil and gas industry. In the face of increasing global energy consumption, as well as the need to reduce environmental stress, enhanced Oil Recovery (EOR) and Carbon Capture and Sequestration (CCS) may be a future part of hydrocarbon production in a more sustainable manner. With EOR, more hydrocarbons can be obtained from a single well by injection of brine, polymer, or surfactant, thereby reducing the need for additional drilling. In CCS, carbon dioxide is re-injected into the reservoir so that the carbon dioxide does not enter the atmosphere. The flow and dispersion of fluids in rock also drives the geothermal extraction process, which is another promising source of renewable energy. Understanding fluid flow properties plays a key role in the fine tuning and deployment of such techniques.

However, many factors affect fluid flow, including, for example, fluid density and viscosity, interfacial tension, fluid flow rate, surface wettability, and pore geometry. Laboratory testing of these factors requires a significant amount of time and is very expensive. Furthermore, the number of samples that can be processed is relatively limited due to the time and expense required to perform each test. There is a need to provide information more quickly in order to make more timely decisions.

Many in the industry are turning to simulations of fluid flow in reservoirs in order to make informed decisions in a timely manner. Performing pore scale flow simulation based on an understanding of the dynamics of fluid displacement at the pore scale may facilitate computation of darcy scale flow parameters used as inputs, for example, in reservoir simulation.

Direct flow simulation is widely used to predict fluid flow and transport through complex porous rocks. One type of simulator is based on the lattice-Boltzmann method (LBM) which simulates collisions and flow of microscopic particles and then evaluates macroscopic pressure gradients and velocities, whereby the permeability of the porous material can be estimated (Chen et al, "lattice Boltzmann method for fluid flow (Lattice Boltzmann method for fluid flows)") Fluid mechanics annual book (Annu_Rev\u) Fluid Mech)》30:329-64; 1998). For example, alpak et al ("predicting fluid topology and relative permeability (Prediction of fluid topology and relative permeability in imbibition in sandstone rock by direct numerical simulation) in sandstone imbibition by direct numerical modeling")Water resource development (Advances_in_Water_Resources)》122:49-59;2018 and "two-phase viscous capillary flow (Direct simulation of pore-scale two-phase visco-capillary flow on large digital rock images using a phase-field lattice Boltzmann method on general-purpose graphics processing units)" on a large digital rock image using phase field lattice boltzmann method directly simulating pore scale on a general graphics processing unitComputing earth science23:849-880;2019 Describes an energy-based LBM (eLBM), which is composed of three componentsTwo-phase flow simulation system consisting of individual modules: the device comprises a forced discharge simulation module, a forced imbibition simulation module and a steady-state relative permeability calculation module.

Before such simulators can be effectively used as predictive tools, their output should be verified against experimental data. Early studies have compared Nuclear Magnetic Resonance (NMR) and Magnetic Resonance Imaging (MRI) measurements of flow in porous materials (flow MRI) with direct flow simulation in porous materials. However, while good qualitative agreement is shown between experimental velocity graphs and simulations, LBM simulations tend to predict large axial (in the surface flow direction) velocities too high, but low axial velocities too low.

MRI and flow MRI are useful tools for characterizing microstructure and mass transport properties in porous media. Such data is typically obtained at a spatial resolution of several hundred microns that is too coarse to adequately capture pore scale information in porous rock. Therefore, there is a need to increase the spatial resolution of MRI acquisitions by 1 to 2 orders of magnitude to explore structural transport relationships at spatial resolutions closer to that of the relevant pore scale. However, 3D MRI acquisition at this spatial resolution is impractical due to the lengthy acquisition time. Rapid 3D MRI acquisition methods, such as RARE (rapid acquisition with relaxation enhancement), alleviate this problem by significantly reducing acquisition time.

De Kort et al ("undersampling and compressed sensing of 3D spatially resolved Displacement propagators in porous media Using APGSTE-RARE MRI (render-sampling and compressed sensing of 3D spatial-resolved displacement propagators in porous mediausingAPGSTE-RARE MRI))Magnetic Co-ordination Vibration Imaging (Mag Res Imaging)56:24-31;2019 A schematic of a flow MRI pulse sequence for acquisition of spatially resolved propagators is presented. The sequence combines an alternating pulse gradient stimulated echo (APGSTE) sequence for encoding displacement with a RARE imaging sequence used in a velocity imaging sequence.

Karlsons et al ("identification of sampling patterns for high resolution compression-sensing MRI of porous materials: from"Learning" in X-ray microcomputerized tomographic data (Identification of sampling patterns for high-resolution compressed sensing MRI of porous materials: "Learning" from X-ray micro-computed tomography data) "Microscope journal of academic journal j.microsc.)276:63-81;2019 A data sampling pattern for implementing an undersampling scheme for compressed sensing MRI (CS-MRI) acquisition is disclosed.

Digital petrophysical is a technique that has been developed to provide faster, more and less expensive analysis of hydrocarbon-bearing formation rock to determine key petrophysical properties of the rock. Digital petrophysical utilizes digital images of formation rock to simulate petrophysical properties of pore dimensions and predict properties of complex rock. The main focus of DR technology is to accurately simulate fluid flow behavior within the pore space in order to predict petrophysical properties of the sample, such as permeability.

There is a need for a method for calibrating direct flow simulation of rock samples using flow MRI images acquired at a spatial resolution compatible with pore-scale direct flow simulation.

Disclosure of Invention

According to one aspect of the present invention there is provided a method for calibrating direct flow simulation of a rock sample, the method comprising the steps of: providing a 3D image of the rock sample, wherein the 3D image is comprised of a plurality of voxels; generating a segmented structure image of the rock sample from the 3D image by selecting each voxel of the plurality of voxels to represent a void space or a solid material; simulating fluid flow on the segmented structure image using direct flow simulation; generating a 3D spatially resolved fluid velocity map of one or more fluid phases at a pore scale resolution using pulsed field gradient magnetic resonance imaging; and comparing the simulated fluid flow and the 3D spatially resolved fluid velocity map; thereby calibrating the direct flow simulation on the rock sample.

Drawings

Several of the figures of the drawings showing non-limiting examples of the invention are performed in the original submitted color and are not easily converted to black and white in a manner that shows the information being conveyed. The color drawings will also be submitted for publication in non-patent locations.

The process of the present invention may be better understood by reference to the following detailed description of preferred embodiments and the accompanying drawings referred to therein, in which:

FIG. 1 is a block flow diagram illustrating the steps of one embodiment of the present invention;

FIG. 2 is a block flow diagram illustrating steps of another embodiment of the present invention;

FIG. 3A is a schematic diagram of an exemplary PGSE-RARE flow MRI pulse sequence that may be used to acquire velocity maps in one embodiment of the invention;

FIG. 3B is a schematic diagram of an exemplary APGSTE-RARE flow MRI pulse sequence that may be used to acquire spatially resolved propagators in another embodiment of the invention;

FIGS. 4A and 4B are 3D sections extracted from microscopic CT images of Ketton rock samples and Estaillades rock samples prepared in a non-limiting embodiment of the invention;

FIG. 5A is a pore scale 3D magnitude flow chart showing the main flow channel through a Ketton rock sample in a non-limiting example of the invention;

FIGS. 5B and 5C are graphical velocity profiles of the z-velocity component and the y-velocity component, respectively;

FIGS. 6A-6C provide visual comparisons of direct flow simulation of Ketton rock samples and color coding of co-registered images in a non-limiting example of the invention, where μCT image density is gray scale (ranging from low intensity black to high intensity white) and speed is color coding scale (ranging from-2 mm/s blue to 4mm/s red);

FIG. 7 is a graphical depiction of the velocity profile of a direct flow simulation and MRI velocity map of a Ketton rock sample in a non-limiting example of the invention;

FIG. 8 is a 2D histogram showing voxel-by-voxel correlation between a direct flow simulation of Ketton rock samples and an MRI velocity map in a non-limiting example of the invention;

FIGS. 9A and 9B are graphical pore scale comparisons of direct flow simulation and MRI velocity maps of Ketton rock samples in a non-limiting example of the invention;

FIG. 10 is a color-coded 2D slice image of a 3D co-registered image of Estalaillades rock sample from a non-limiting example of the present invention, where the μCT image density is gray scale (ranging from low intensity black to high intensity white) and the MR signal intensity is color-coded scale (ranging from low intensity blue to high intensity green);

11A-11I provide visual comparisons of LBM simulations of Estaillades rock samples in non-limiting examples of the present invention and color coding of MRI velocity maps, where μCT image density is gray scale (ranging from low intensity black to high intensity white) and velocity is color coding scale (ranging from-2 mm/s blue to 5mm/s red);

FIGS. 12A and 12B are graphical pore scale comparisons of direct flow simulation and MRI velocity maps of Estaillades rock samples in a non-limiting example of the present invention; and is also provided with

FIG. 13 is a graphical comparison of the flow rates of each image slice in the middle region of an Estaillades rock sample in a non-limiting example of the present invention.

Detailed Description

The present invention provides a method for quantitatively visualizing fluid transport properties of rock. According to the invention, a direct flow simulation of fluid flow is compared to a 3D spatially resolved fluid velocity map to calibrate the direct flow simulation on the rock sample. By better understanding the flow and distribution of fluids over the rock sample, better decisions can be made regarding recovery of hydrocarbons from the subsurface formation, carbon capture and sequestration in the subsurface formation, and/or geothermal extraction processes.

Flow MRI provides the ability to quantitatively and non-invasively measure 3D spatially resolved fluid flow fields in rock samples. By looking at the spatially resolved flow field, it can be determined whether the microscopic flow phenomenon predicted based on core-level petrophysical properties was properly simulated. Furthermore, spatially resolving the flow properties provides a basic insight into the degree of spatial non-uniformity of the flow over the rock sample.

Referring now to fig. 1, one embodiment of the method 10 of the present invention involves providing a 3D image of a rock sample at step 14 and generating a segmented structure image of the rock sample at step 16. The segmented structure image identifies the void space and solid material in each voxel of the segmented structure image. At step 18, fluid flow is simulated on the segmented structure image using direct flow simulation.

Petrophysical properties of rock samples have been measured using digital rock techniques using segmented structure images. While these macroscopic or global properties are critical to any petrophysical analysis, they provide limited insight into the degree of spatial variation of these properties. The non-uniformity of petrophysical properties is a key factor in the efficiency of the hydrocarbon recovery process. The inventors have surprisingly found that digital rock technology is enhanced by spatially resolving fluid flow in a rock sample. Thus, in accordance with the present invention, a 3D spatially resolved fluid velocity map of one or more fluid phases is generated at step 22.

The simulated fluid flow and velocity maps are compared in step 24 to correlate the pore structure of the rock sample with the fluid flow in the pores. By comparing the pore structure to a direct simulation of fluid flow and a 3D spatially resolved velocity map, the direct flow simulation can be calibrated. The comparing step may be performed by visual comparison between the pore structure image and the fluid flow simulation image and the 3D spatially resolved velocity map. In a preferred embodiment, the comparing includes co-registering the simulated fluid flow and the 3D spatially resolved fluid velocity map to produce a co-registered 3D image. The co-registered 3D images are used to calibrate the direct flow simulation at all locations within the rock sample.

In the embodiment of fig. 2, the segmented structure image generated in step 18 is used to generate a 3D spatially resolved fluid velocity map at step 22. In this embodiment, the segmented structure image is used to optimize the data undersampling pattern of the experiment. Compressive sensing is then used to reconstruct undersampled data acquired according to the pattern.

3D image

The 3D image of the rock is obtained from the rock from a formation of which the fluid transport properties are of interest. By way of example, the rock may be sandstone, carbonate, shale, and combinations thereof from a hydrocarbon containing formation. Alternatively, the rock may be from a subterranean formation in which carbon sequestration is considered. The rock may be obtained by conventional means for obtaining rock samples from the formation. In a preferred embodiment, a core sample of rock is obtained by coring a portion of the formation (e.g., the entire core or a sidewall core) from a well within the formation. Alternatively, the rock sample may be obtained from drill cuttings, preferably undisturbed drill cuttings, produced while drilling a borehole in the formation. Rock may be obtained from the same borehole as the electrical property measurements. Alternatively, the rock may be obtained from another wellbore in the same field as the wellbore that produced the electrical property measurement.

The rock sample should be of sufficient size to obtain a 3D image of sufficient volume at the scale at which the image is generated. In particular, the rock sample should be of sufficient size such that at the scale or field of view of the image to be generated, the properties of the bulk of the sample are superior to those of the edges of the sample.

A 3D image consisting of a plurality of voxels is obtained from a rock sample. A 3D image of the rock may be obtained using pore scale imaging techniques. The 3D image of the rock may be obtained by, for example, but not limited to, scanning Electron Microscopy (SEM), X-ray computed tomography, acoustic microscopy, magnetic resonance imaging, and the like. X-ray computed tomography includes, but is not limited to, X-ray microcomputerized tomography (micro CT) and X-ray nano computerized tomography (nano CT). Most preferably, a 3D image of the rock is obtained by micro-CT to provide a sufficient field of view of the rock to avoid edge porosity distorting the resulting image, as well as reducing the scan time and computational requirements that would be required for higher resolution tomography (e.g., nano-CT).

A 3D image of rock obtained by aperture scale imaging techniques is made up of a plurality of voxels, wherein the volume defined by each voxel represents the maximum resolution of the image. The resolution of the image should be selected to provide a voxel size at which the primary pore throat for fluid flow in the rock is sufficiently resolved and at which a sufficient field of view is provided to represent the entire rock of the fluid transport characteristic to be analyzed.

The resolution of the micro CT image may be selected based on the size of the rock sample, the relative average pore size of the rock type, the time required for imaging, and the computational power required to store and perform further computational activities on the image data. The range of pore-scale resolution of the micro-CT image may be 0.1 μm per voxel ³ To 30 μm ³ . For sandstones, the micro-CT image is preferably in the range of 1 μm per voxel ³ To 25 μm ³ More preferably 2.5 μm per voxel ³ To 15 μm ³ Is generated. For carbonate rock, the resolution of the micro CT image is preferably in the range of 0.5 μm per voxel ³ To 20 μm ³ More preferably 1 μm per voxel ³ To 10 μm ³ Is generated. For shale, the resolution of the micro CT (or nano CT) image is preferably in the range of 0.1 μm per voxel ³ To 10 μm ³ More preferably 0.5 μm per voxel ³ To 5 μm ³ Is generated.

In a preferred embodiment, the acquired image may be processed to reduce noise and image artifacts. Noise may be filtered from the acquired image by filtering using a local mean filter to reduce noise. The imaging artifacts dominant at the outer edges of the acquired image can be reduced by processing the image while excluding the outer edges of the image.

Segmentation

The 3D image obtained for the rock is processed to segment voxels of the image into voxels representing either a void space in the rock or a solid material in the rock, thereby producing a binary image in which the void voxels have a value of zero and the solid material voxels have a value of one (or vice versa). The image may be a gray scale image and processing voxels of the image to divide the image into voxels representing void space or solid material may be performed by assigning voxels as holes based on a threshold valueA void space or as a name for a solid material, wherein voxels with an image intensity above a threshold may be assigned a value representing a void (or solid material) and voxels with an image intensity below a threshold may be assigned a value representing a solid material (or void). Can be used in Otsu ("threshold selection method from gray level Histogram (A Threshold Selection Method from Gray-level Histogram)")IEEE System, human and control theory journal (IEEE Trans.SMC)》9:62-66;1979 Otsu method described in the art or other threshold calculation algorithm known in the art).

Segmentation algorithms are known to those skilled in the art. Preferably, the segmentation method is selected to identify pore space from the solid matrix. Examples of segmentation methods are described in Otsu ("threshold selection method from gray level Histogram (A Threshold Selection Method from Gray-level Histogram)") IEEE System, human and control theory journal (IEEE Trans.SMC)》9:62-66;1979 A) is provided; andra et al ("digital petrophysical benchmark-Part II: computationally efficient Properties (Digital Rock Physics Benchmarks-Part II: computing Effective Properties)")Meter (Meter) Computer and geology science (Computers and Geosciences)50:33-43;2013 A) is provided; the effects of Saxena et al ("image segmentation and voxel size on the efficient transport and elastic properties of the micro CT calculation (Effect of Image Segmentation)&Voxel Size on Micro-CT Computed Effective Transport&Elastic Properties)″Ocean and Petroleum geology (Marine and Petroleum Geology)86:972-990;2019 A) is provided; and Chuang et al ("Fuzzy C-means clustering (Fuzzy C-Means Clustering with Spatial Information for Image S egmentation) using spatial information for image segmentation)"Medical imaging of the book And graphics (Comput. Med. Imaging graph.)30:9-15;2006 Described in (c). Those skilled in the art will understand the desired choice of segmentation. Segmentation using a segmentation algorithm is preferably automated using a data processing system.

In a preferred embodiment, the 3D image is segmented by a watershed-based segmentation algorithm (Beucher et al, "morphological method for segmentation: watershed transformation (The morphological approach to segmentation: the watershed transformation)", E.R. Dougherty (editorial), "mathematical morphology image method (Math. Morph. Image Process.), new York Marcel Dekker Inc., new York, 1993: pages 433 to 481).

Direct flow simulation

According to the invention, the fluid flow is simulated on the segmented structure image using direct flow simulation. Direct flow simulation includes, for example, but is not limited to, finite difference methods, finite element methods, finite volume methods, and lattice boltzmann methods. In a preferred embodiment, the fluid flow is simulated with an LBM simulator. Examples of LBM simulators include, but are not limited to, energy-based LBM (eLBM) simulators and Multiple Relaxation Time (MRT) LBM simulators.

Spatially resolved fluid velocity mapping

The fluid velocity profile is obtained experimentally. According to one embodiment of the invention, the segmented structure image is further used to generate a 3D spatially resolved fluid velocity map of one or more fluid phases at a pore scale resolution using pulsed field gradient magnetic resonance imaging.

Preferably, the 3D spatially resolved fluid velocity map is generated by: a spatially resolved fluid flow encoded phase difference map of each of the fluid phases in the fluid phase in a flow state is reconstructed using Compressive Sensing (CS), and a flow velocity map of each of the fluid phases is calculated. The flow rate map is preferably calculated by: correction is applied to any spurious phase angle shifts introduced by the relaxation enhanced Rapid Acquisition (RARE) MRI pulse sequence. Preferably, the correction is applied by subtracting the static phase difference from the phase difference measured in the flowing state.

MRI has been widely used to study porous materials, including porous core plugs, because it provides a direct and non-invasive way to obtain quantitative chemical, microstructure, and transport related information within fluid saturated optically opaque materials. Due to the relatively low sensitivity of Magnetic Resonance (MR) methods, such MRI data are typically acquired with a spatial resolution on the order of hundreds of microns. However, the pore size in sedimentary rock is typically less than 100 μm. There is a need to increase the spatial resolution of MRI measurements in order to be able to obtain local structural transport correlations with a spatial resolution similar to that of the relevant pore scale.

CS in MRI is based on requirements such as (1) aliasing artifacts (e.g., sample data) in linear reconstruction must be incoherent and noise-like; (2) the desired image exhibits transform sparsity; and (3) reconstructing the image using a non-linear algorithm that implements sparsity and consistency with the acquired k-space data.

As discussed in Karlsons et al (2019), 3D MRI images can be acquired using a combination of fast and undersampled MRI data acquisition and CS at a spatial resolution of 17.6 μm (carbonate) that is at least an order of magnitude higher than the resolution of conventional MRI and is a resolution that can clearly discern pore scale features in many rocks. CS significantly reduces experimental acquisition time by reconstructing images based on undersampled data space using a priori knowledge.

Investigation of a sparsity-promoting variational approach in CS flow MRI (Benning et al "phase reconstruction from velocity-encoded MRI measurements" (Phase reconstruction from velocity-encoded MRI measurements-A survey of sparsity-promoting variational approaches) "Magnetotonresonance journal (J Magnetic) Resonance)》,234:26-43;2014; lusting et al "sparse MRI: application of compressive sensing for fast MR imaging (spark MRI: the application of compressed sensing for rapid MR imaging) "Medical science MagnetotonMagnetic magazine (Magnetic) Resonance in Medicine)》58:1182-1195;2007 The goal is to recover the phase m of the complex-valued image from the acquired undersampled k-space dataset y. CS enables m, m to be implemented using a variational regularization method that incorporates a priori knowledge about m into the reconstruction process _CS Is a solution to (a). In the case of CS reconstruction, a priori knowledge is that the image can be sparsely represented in the appropriate transform domain.

The type J of regularization function used to map the image into the transform domain depends on the nature of the image to be mapped. For example, a non-smooth regularization such as Total Variation (TV) may be better suited for images with sharp edges, while a smooth regularization such as Daubechies wavelet transform is well suited for images where the pixel intensity variation is more gradual.

In one embodiment, TV regularization based on m finite difference transforms is used for CS reconstruction. Calculating m based on this a priori knowledge using a nonlinear recovery algorithm _CS Is expressed as

Where α is the equilibrium fidelity termRegularization parameters for weights between regularization term αj (m). The value of the parameter alpha is selected based on the Morozov deviation principle, which is written as

Where n is the number of k-space samples. The maximum value of α that satisfies the condition in equation (2) is selected for reconstruction.

Imaging may be performed to selectively measure various fluids in a core sample, such as hydrocarbons (e.g., crude oil or dodecane) and aqueous fluids (e.g., water, brine, etc.). Such techniques may be used to image various fluids in the formation, alone or in combination. In particular, imaging may be used to distinguish aqueous fluids from hydrocarbons in a core sample. These images may be used, for example, to characterize fluid parameters such as flow rate and type of hydrocarbon produced. Information collected from such imaging may be used, for example, to identify a particular fluid, to image the fluid alone, to evaluate formations containing the fluid, to determine downhole parameters, to detect valuable hydrocarbons, to provide information for planning oilfield operations, and so forth. Flow imaging is performed while one or more fluids are pumped through the rock sample, alone or in combination. The static state is imaged by interrupting the fluid pumping.

Flow MRI is an integration of pulsed field gradient nuclear magnetic resonance (PFG NMR) for quantifying fluid displacements with MRI that spatially resolves these displacements. Preferably, PFG NMR is combined with a RARE MRI pulse sequence.

A schematic illustration of a velocity encoded flow MRI pulse sequence for acquiring velocity maps is shown in fig. 3A, which will be discussed in more detail in an embodiment. The pulse sequence combines a PGSE sequence for encoding velocity with a RARE imaging sequence for spatially encoding the measured velocity. The RARE imaging sequence is well suited for CS applications and imaging porous rock that is fluid saturated. PGSE sequences for velocity encoding are preferred because the observation time is short compared to the signal decay rate in the transverse plane. Thus, the signal loss is small compared to using stimulated echo sequences which inherently suffer from 50% signal loss. In a preferred embodiment of the invention, both RARE imaging sequences and CS methods are utilized to accelerate acquisition and enable a 35 μm spatial resolution 3D velocity map to be obtained.

In another embodiment of the invention, spatially resolved displacement propagators are acquired. The displacement propagator is a probability distribution of molecular displacement Wherein->Is observed at an observation time delta (typically in a few ms up to the spin-lattice relaxation time T of the fluid ₁ In the range of (a), the spin-lattice relaxation time can be as high as a few seconds) the probability of the spin moving over the dynamic displacement r=r' -R. By spatially resolving the propagators, information about local time-dependent flow dispersion due to molecular self-diffusion of fluid molecules across the streamline, as well as measurements of flow velocity, can be obtained.

A schematic diagram of the acquisition of a spatially resolved train of pulses of propagators is shown in fig. 3B, which will be discussed in more detail in an embodiment. The sequence will be used for coding the displacementAn alternate pulse gradient stimulated echo (APGSTE) sequence is combined with a RARE imaging sequence used in a velocity imaging sequence. The stimulated echo sequence provides a longer observation time for the propagator measurements. APGSTE or Cotts 13-interval based on stimulated echo sequence (Cotts et al "pulsed field gradient stimulated echo method (Pulsed field gradient stimulated echo methods for improved NMR diffusion measurements in heterogeneous systems)" for improving NMR diffusion measurements in non-uniform systemsMagnetotonresonance magazine (J Magnetic Resonance)》83:252-266;1999 The sequence utilizes alternating pulse gradients to reduce signal loss due to motion in the background field gradient and is well suited for measuring flow in porous media. Alternatively, a Pulse Gradient Stimulated Echo (PGSE) sequence may be used to encode the displacement.

In one embodiment of the invention, a single average velocity is spatially resolved for each of the voxels within the image-a method known as velocity mapping. One method that may be used for acquisition velocity mapping is a Pulse Gradient Spin Echo (PGSE) sequence. Preferably, PGSE is used with the RARE imaging sequence.

Undersampled k-space datasets in the flow state and in the static state are reconstructed separately to produce complex-valued MRI images. From these phase diagrams, a velocity map is then calculated as described above. Preferably, the background noise is reduced, more preferably zero, by generating a binary mask from the static MRI dataset and multiplying the mobile phase map by the binary mask.

A 3D velocity map, and preferably a static velocity map, is generated from the mobile phase to show the main flow path through the rock sample. The 3D velocity map provides a visualization of the non-uniformity of the rock sample and the presence of higher velocity flow channels. Further, areas of rock having low fluid mobility and/or retained fluid may be visualized. Finally, the negative velocity accounts for the occurrence of backflow through porous rock and is believed to be caused by the recirculation flow pattern or vortex-like structure at the junction of streamlines immediately adjacent the surface.

In accordance with the method of the present invention, a high-resolution MRI velocity map capable of non-invasively capturing flow information in the pore space is compared to a direct flow simulation run on a segmented 3D muct image plug.

In a preferred embodiment of the method of the present invention, the results of the flow simulation are compared to the MRI flow map by co-registering and coarsening the simulation down to the same resolution as the flow MRI data. Co-registering the simulated fluid flow with the 3D velocity map allows the simulated fluid flow to be calibrated with the experimentally generated velocity map. This is particularly advantageous for rocks with structural non-uniformities and/or fluid flow in the pores below the resolution of the 3D image.

X-ray muct images of small rock bolts are acquired at a resolution of a few microns (i.e., an order of magnitude higher than the resolution currently available for flow MRI). However, especially for heterogeneous carbonates, the smallest (micro) pores are still significantly smaller than the currently available voxel size of X-ray μct. The region of rock where microporosity is present will have a lower average X-ray absorption because the overall lower mass density of that region is such that such microporous regions cannot be distinguished from non-porous regions having a lower mineral density so that the X-ray absorption is effectively the same, unless the rock is saturated with doping fluid.

The challenge of conventional methods is to understand the extent to which these smaller pores actually contribute to the overall flow of fluid through the rock, and thus the accuracy of the simulation of the flow field of the pore space derived based on μct. Although in many cases it is not desirable that the microporosity of the rock carries significant flow, some regions of microporous carbonate rock may be connected by microporosity alone, in which case the contribution of microporosity to flow may become significant.

According to the invention, the comparison of the 3D velocity map with the simulated fluid flow makes it possible to distinguish between the two non-invasively. MRI is only sensitive to fluids and will therefore show only fluid saturated microscopic pores, not solid rock areas where no fluid is located. In an embodiment of the invention, wherein the comparing step comprises co-registering the 3D velocity map with the simulated fluid flow, the co-registered 3D image providing an average fluid flow velocity for each voxel in the segmented structure image.

Images depicting fluid flow through the rock are co-registered in a predetermined flow direction to determine a fluid flow distribution over the rock sample. Image co-registration may be performed using software tools including, for example and without limitation, aviZO ^TM 2019.4 The "Register Images" module in (FEI visualization sciences, USA, FEI Visualization Sciences Group) USA) is completed.

In one embodiment of the invention, the segmented structure image is co-registered with a spatially resolved fluid velocity map. In another embodiment of the invention, a rigid transformation is used to optimize alignment of the co-registered structure and velocity map. The transformed image is then co-registered with the direct flow simulation.

In yet another embodiment of the present invention, lanczos filters are used (Duchon, "Lanczos filtering in one and two dimensions (LanczosFiltering in One and Two Dimensions)")Using meteorological and climatic impurities Shi (j. Appl. Metaorol.)18:1016-1022;1979 Resampling the registered images onto the MRI data coordinate system, which enables the MRI flowgram and 3D images (including LBM simulations) to be visualized and compared on the same grid.

Examples

The following non-limiting examples of embodiments of the methods of the invention claimed herein are provided for illustrative purposes only.

The rock samples used in this example were Ketton limestone core plugs (3.84.+ -. 0.01mm in diameter and 11.10.+ -. 0.37mm in length) and Estaillades limestone core plugs (3.92.+ -. 0.02mm in diameter and 12.65.+ -. 0.13mm in length). The samples were dried in an oven at 70 ℃ for 12h.

Example 1 structural imaging

Using Bruker SKYSCAN 1172 ^TM A Micro CT scanner (Bruker Micro-CT, belgium) images the dried Ketton sample at an isotropic resolution of 5.00 μm. Performed using a 60kV source voltage, a 165 μA source current, and an Al (0.5 mm) filterAnd (5) collecting row images. 802 projection images were acquired by rotating the sample within 200.5 ° at 0.25 ° angular increments, with 10 scans per angular increment, resulting in a total acquisition time of 11.4 h. Using NRECON ^TM Software (Bruker, v1.6.8.0) reconstructs the projection images to obtain 2666 cross-sectional slices.

Using Bruker SKYSCAN 2214 ^TM A Micro CT scanner (Bruker Micro-CT, belgium) images the dried Estaillades sample at an isotropic resolution of 3.00 μm. Image acquisition was performed using a source voltage of 90kV, a source current of 70 ua and an Al (1 mm) filter. To image the entire sample, acquisitions are performed at 5 different scan positions along the longest dimension of the sample. For each position 3601 projections were obtained by rotating the sample in angular increments of 0.1 ° over 360 °, with 6 scans per angular increment yielding an acquisition time of 4.26 h. Thus, the total acquisition time was 21.3h. Stitching together the projected images from all 5 scan positions and using NRECON ^TM Software (Bruker, v1.6.8.0) was reconstructed to obtain 4319 cross-sectional slices.

A 3D CT image of the rock sample is generated by stacking all 2D cross-sectional slices consecutively. For the Ketton rock sample and the Estaillades rock sample, 3D slices 20 extracted from the μCT images of the rock samples are shown in FIGS. 4A and 4B, respectively.

Then in AVIZO ^TM 2019.4 The 3D muct image was segmented (FEI visualization sciences, usa) using a watershed-based segmentation algorithm (Beucher et al, 1993). A binarized structural image is generated by assigning all non-0 phases to 1 phase, which is intended to represent an impermeable solid. The pore space of the binarized structural image was separated into individual pores using a Chamfer distance transform and a marker-based watershed algorithm.

The rock sample was vacuum saturated with deionized water. Pore network properties were extracted from the μct images of the Ketton core plug and the estailades core plug as shown in table 1. Properties are extracted from the connected pore space images. The spatial resolution of the image was 7 μm for the Ketton sample and 3 μm for the Estaillades sample.

TABLE 1

Example 2 fluid flow imaging

For flow MRI experiments, samples were placed in an Adtech Fluorinated Ethylene Propylene (FEP) heat shrink tubing to connect the samples to the inlet and outlet of the FEP flow line and provide constraints for creating a surface flow direction (z) that represents the net total flow. Although fluid may also flow in the x-plane and y-plane, the net total flow in the x-direction and y-direction will be 0, as fluid cannot flow out the sides due to restriction. In the case of the Estaillades sample, two layers of Teflon tape 75 μm thick were applied around the plug prior to placing the sample in the heat shrink tubing to minimize the diversion of fluid through the grooves and cracks on the surface of the rock during the flow experiments. Using Vindium VP-6 ^TM The metering pump applies a constant water flow rate.

Three pore-scale MRI datasets were acquired: a velocity map and two spatially resolved propagators, which have different observation times (delta), as described in more detail below. In Bruker BIOSPIN AVANCE III HD ^TM Experiments were performed on a spectrometer controlled 7.0T vertical bore magnet. Using a magnetic field having 2.9T m in three orthogonal x, y and z directions ^-1 Bruker MICRO5 for maximum gradient intensity of (2) ^TM A three-axis gradient system to achieve spatial resolution. Tuning to 299.84MHz ¹ H) A 8mm radio frequency (r.f.) saddle coil for spin excitation and signal detection.

Example 2A-MRI dataset 1: velocity map

Undersampled velocity maps were acquired using PGSE-RARE experiment 32, schematically shown in fig. 3A. PGSE-re experiments 32 combine PGSE sequences 34 with re imaging sequences 36 for use in velocity and spatial encoding, respectively. The PGSE 34 velocity encoding is achieved by a pair of intensity unipolar gradients g durations delta separated by an observation time delta. Frequency of useCoding gradient G _z Phase encoding gradient G _x And G _y To implement the RARE 36 spatial coding. RARE factor (N) _RF ) Is the number of 180 refocusing pulses in the echo train. Is the echo interval.

The sampling pattern for the acquisition velocity map is generated using a μCT-based variable density sampling (μCT-VDS) strategy described in Karlsons et al (2019). In this embodiment, for Ketton samples, the method is used to generate an optimized k-space sampling pattern with a k-space sampling fraction of 0.25. The duration of the hard 90 ° excitation and 180 ° refocusing r.f. pulses were 6.0 μs and 12.0 μs, respectively. Use as N _RF RARE factor of=8 (N _RF ) And is t _e Echo interval of =2.2 ms (t _e ) To acquire N of k-space for each excitation pulse _RF A continuous line. To achieve speed encoding, the amplitude g=2.0T m is used ^-1 Is a pair of unipolar gradients g (g) _i ＝4.0T m ^-1 ) And duration δ=0.132 ms, the pair of unipolar gradients and duration being separated by an observation time Δ=4 ms. For averaging and recirculating the delay (t _RD ) Is L _RD 32 scans of =1.1s, the acquisition time of the velocity profile was 20h. The velocity is encoded along the z-direction, i.e. along the surface flow direction. Images in a static state are also acquired to correct for speed shifts. Thus, the total acquisition time of the velocity map with the velocities encoded in the surface flow direction is 40h. Images were acquired with a field of view (FOV) of 13.5mm x 4.5mm and 384 voxels x 128 voxels in the frequency encoding (z) and phase encoding (x, y) directions, respectively, yielding a 3D velocity map with an isotropic resolution of 35.2 μm.

For Estaillades samples, undersampled velocity maps were acquired using the same PGSE-RARE pulse sequence. A variable density sampling based on μct with a k-space sampling fraction of 0.3125 is used to generate the undersampling scheme. Hard 90 ° excitation and 180 ° refocusing r.f. pulses were used for 6.6 μs and 13.2 μs, respectively. The number of echoes in the RARE echo chain is NR _F Echo interval t=8 _e =2.2 ms. Speed imaging parameter g _i ＝2.3T m ^-1 Delta=0.132 ms and delta=4 ms. To be used fort _RD 32 scans were acquired for signal averaging, giving an acquisition time of 31.3h =1.375 s. The total acquisition time of the offset corrected velocity map is about 63h. A velocity map was acquired with a FOV of 13.5mm x 4.5mm and 384 voxels x 128 voxels in z direction, x direction and y direction, yielding a 3D velocity map with an isotropic resolution of 35.2 μm.

Example 2B-MRI datasets 2 and 3: propagator

Spatially resolved propagators were acquired using an APGSTE-RARE experiment 42 schematically shown in fig. 3B. The APGSTE-RARE experiment 42 consists of a 13-spaced APGSTE sequence 44 and a RARE imaging sequence 36. Displacement encoding is achieved by a gradient g, which can be applied in any of three orthogonal directions (z, x or y). A displacement encoding gradient is applied parallel to the surface flow direction (z). Encoding gradient G by frequency _z And phase encoding gradients Gx, y to achieve spatial encoding.

APGSTE-RARE experiment 42 was performed with the following parameter settings: the duration of the hard 90 ° excitation and 180 ° refocusing r.f. pulses were 6 μs and 12 μs, respectively. Number of echoes N acquired in each echo train _RF =8, and t _e ＝4.2ms。

Two different spatially resolved propagators are acquired at different observation times Δ, i.e., Δ=150 ms and Δ=900 ms. Applying a gap flow velocity v _i Scaling as reciprocal of the observation time such that the total gap displacement v _i Delta is the same for both experiments and self-diffusion effects can be observed. For an observation time of 150ms, 0.03ml min was applied ^-1 (v _i =91±3ft days ^-1 (0.32±0.01mm s ^-1 ) At a fluid flow rate of 0.005ml min ^-1 (v _i =15±1ft days ^-1 (0.05±0.003mm s ^-1 ) An experiment with an observation time of 900ms was performed at a flow rate.

The displacement encoding is by a duration of delta/2=0.14 ms and maximum amplitude g=2.8T m for delta=150 ms ^-1 And g=1.3T m for a maximum amplitude of Δ=900 ms ^-1 Is achieved by a two pair bipolar gradient g.

3D spatially resolved propagators were acquired with a FOV of 13.5mm x 4.5mm and 144 voxels x 48 voxels in the frequency direction (z) and two phase encoding directions (x and y), respectively, yielding a 3D image with an isotropic resolution of 93.8 μm. In q-space, 65 points are acquired, resulting in a FOV in the displacement direction of 1.0mm. The displacement direction is selected to be parallel to the surface flow direction (z).

For 3D spatially resolved propagator experiments, 3D q-k is required _{Phase (C)} Spatial undersampling pattern. Thus, the shape of the probability density function (pdf) describing the undersampled mode is based on (1) for 2Dk _{Phase (C)} The μCT-VDS method of the space and (2) the shape of the 1D q space derived from the spatially unresolved propagation sub-experiment. The two pdfs are multiplied to obtain a 3D probability distribution function based on which a sampling pattern is generated with a sampling score of 0.25.

Using t _RD The mean value of 16 signals was acquired with a sampling fraction of 0.25 in k, q space for 2s, the acquisition time of spatially resolved propagators was 4 days for the Δ=150 ms experiment and 5 days for the Δ=900 ms experiment. This is 32 times faster than a fully sampled simple spin echo MRI experiment given the signal-to-noise ratio.

For both the velocity map and the spatially resolved propagation sub-experiment, the sample remains in the same position.

Example 3-image processing

Phase reconstruction of MRI measurements from velocity encoding using an additional MATLAB toolbox, specifically object-oriented math for inverse problem (OOMFIP) (Benning et al "investigation of sparsity-promoting variational methods (Phase reconstruction from velocity-encoded MRI measurements-A survey of sparsity-promoting variational approaches)") Journal of magnetic resonance (J.Magn.Reson.)238:26-43;2014 CS reconstruction).

Separate undersampled k-space datasets in the flow state and in the static state are separately reconstructed to produce complex-valued MRI images. A velocity map is then calculated from these phase maps using standard procedures. The phase maps of the odd and even echoes are calculated separately and the resulting phase difference maps of the odd and even echoes are recombined by voxel-wise averaging to improve the accuracy of the velocity map. The amplitude images of the static MRI dataset are then averaged to generate a binary mask, which is multiplied by the velocity map to zero background noise in the velocity map.

Fig. 5A is a pore size 3D magnitude flow chart showing the main flow path through a Ketton rock sample. The graph shows the flow rate v applied from the Ketton sample _i =91±3ft days ^-1 (0.32±0.01mm s ^-1 ) An x-velocity component, a y-velocity component, and a z-velocity component of the water stream. The velocity profiles of the z velocity component and the y velocity component are shown graphically in fig. 5B and 5C, respectively. The x-velocity component has a similar velocity profile as the y-velocity component.

Two characteristics were observed by visual inspection of the flow chart of fig. 5A. First, it can be seen that most of the water in the pore space of the rock sample has a typical velocity < 0.2mm s ^-1 Is not limited, and has a low mobility. Approximately 70% of voxels in this region of the rock contain velocities below about 0.5mm s ^-1 Is a fluid flow of average absolute flow velocity. Second, the flow is non-uniform and concentrated in several high-velocity flow channels.

As determined from fig. 5B and 5C, v _z The distribution of (c) exhibits a long positive tail, at v, indicative of stagnant water _z ＝0mm s ^-1 Significant intensity at the point, and a substantial amount of negative flow or backflow, as represented by the negative velocity. The long positive tail illustrates the non-uniformity of flow through the rock sample, with the highest velocity extending to v _z About 0.08mm s ^-1 About 50 times the modal speed at that time. The average velocity in the surface flow direction (z) is calculated as v _z ＝0.32mm s ^-1 (91 ft day) ^-1 ) This is very consistent with the gap flow rate applied in the macropores. v _y Is distributed in v _y ＝0mm s ^-1 (i.e., zero net flow) distributed symmetrically around the circumference at speeds up to about + -2 mm s ^-1 V of (2) _y . Discovery v _y The positive and negative parts of the distribution of (c) are well described by a single exponential function.

EXAMPLE 4 direct flow simulation

Two LBM algorithms were used to simulate single phase flow in rock. For the Ketton samples, the fluid flow was calculated using an energy-based LBM (eLBM) simulator that utilized a Multiple Relaxation Time (MRT) model as the momentum balance solver. The ehbm is operated in single-phase mode by initializing the pore space with a single-phase fluid and injecting the same phase at the inlet. Therefore, a diffusion interface does not appear in the simulation. In other words, eLBM can be considered as MRT-LBM with constant velocity entry and no gradient exit boundary conditions.

For Estaillades samples, the flow was simulated using the MRT-LBM algorithm based on the same MRT solver.

For Ketton, images of size 1671 voxels x 643 voxels (in z-direction, x-direction, and y-direction) were simulated at an isotropic spatial resolution of 7.0 μm, and for estailades images of size 4279 voxels x 1500 voxels (in z-direction, x-direction, and y-direction) were simulated at an isotropic spatial resolution of 3.0 μm.

To be able to quantitatively compare LBM and MRI speed data, the speed of the high resolution LBM simulation is scaled such that the average slice-by-slice flow rate of the LBM simulation is equal to the applied flow rate of 0.03ml min ^-1 . The flow rate for each slice is calculated by multiplying the sum of all velocities within the slice by the area of each voxel. In the case of Estaillades, the high-resolution analog data is first coarsened (Lanczos filter) to a spatial resolution of 6.0 μm for more efficient image registration and data manipulation. For all quantitative analyses, the analog data was downsampled to 35 μm spatial resolution (Lanczos filter) to enable direct comparison with the acquired MRI speed data.

EXAMPLE 5 Co-registration

Use of "Register Images" modules in AVIZO ^TM Direct flow simulation images and velocity maps of the flow through the rock acquired in the surface flow direction (z) are co-registered.

The purpose of image co-registration is to spatially align the MRI velocity map with the μct-based LBM simulation so that the flow fields of the two methods can be compared qualitatively and quantitatively. An MRI intensity image obtained by calculating the magnitude of the complex-valued static MRI dataset and a μct gray intensity image used to generate the pore space for simulation are used as inputs for image registration.

First, a μCT intensity image (model image) is manually pre-aligned with an MRI intensity image (reference image). Next, a rigid transformation with standardized Mutual Information is used (Plurim et al, "Mutual Information based medical image registration: investigation (Mutual-Information-Based Registration of Medical Images: A surveyy)")IEEE medical image journal (IEEE T) Med Imaging)》22:986-1004;2003 and Studholme et al, "overlap invariant entropy measure of 3D medical image alignment (An overlap invariant entropy measure of 3D medical image alignment)"Pattern recognition (Pattern Recognit)32:71-86;1999 As a measure of the goodness of image alignment, to optimize the alignment of the model image with respect to the reference image. In the optimization procedure, the data is downsampled and image registration is performed step by step with increasing spatial resolution to achieve more efficient image registration. The broad direction and quasi-newton optimizers are used for co-registration of coarse and finer resolution images, respectively. After the image co-registration process is completed, the resulting transformation of the μCT intensity image is applied to LBM simulation.

The co-registered μct intensity image and LBM simulation are resampled onto the MRI data coordinate system using Lanczos filters to enable the MRI flowgram and μct-based image (including LBM simulation) to be visualized and compared on the same grid.

Fig. 6A-6C provide visual comparisons of 2D (xy) image slices extracted from LBM simulations with co-registered images of keton rock samples. Fig. 6A is extracted from a 3D high resolution LBM simulation at a spatial resolution of 7 μm, while fig. 6B is extracted from a co-registered 3D image at a spatial resolution of 35 μm. The images shown in fig. 6A to 6C represent the velocity component in the surface flow direction (z). Good agreement between LBM simulation and MRI speed data was observed. Stored in the MRI flowgram of KettonSubstantially all of the characteristics of the flow behavior in the process are reproduced by LBM simulation. This includes the location of the high-velocity flow channels and the occurrence and location of regions with backflow (i.e., negative (z) velocity). As is evident from the difference plot in fig. 6C, the magnitude of the simulated velocity also agrees well with the acquired MRI velocity, but the LBM appears to underestimate the magnitude of the velocity in the high-velocity flow path of the rock. This observation was consistent with Manz et al ("flow and dispersion in porous media: boltzmann and NMR study of the Lattice (Flow and dispersion in porous media: lattice-Boltzmann and NMR studies)") American chemical industry Society of Engineers (AIChE) J.)》45:1845-1854;1999 Contradiction of reported results, which reports that LBM too high predicts large speeds, possibly due to different LBM codes used to simulate flow. The standard deviation calculated from the middle part of the 3D difference flow graph is 0.37mm s ^-1 . In summary, both the simulation data and the MRI data showed that the flow in the keton rock was highly non-uniform, with the flow carried by several high-velocity flow channels.

To further investigate the consistency between experimental flow data and simulated flow data, a more quantitative analysis was performed by comparing the velocity profile of the two data sets and the voxel-by-voxel velocity. V shown in FIG. 7 _z Is significantly skewed and has a long positive tail extending far beyond the modal velocity of each profile, and a significant degree of reflux, represented by negative velocity. Overall, values greater than v are well predicted by LBM _z ≈0.1mm s ^-1 Is a positive speed of (c). Experimental and simulated distributions also showed that a significant amount of fluid in the interstitial space of the Ketton was stagnant or nearly stagnant, at a velocity of about v _z ＝0mm s ^-1 . This observation is consistent with the visual inspection of the velocity maps in fig. 6A-6C, where the fluid in most of the pores appears to be near stagnant.

However, there is also a significant difference between the two velocity profiles. Two major differences in the shape of the distribution are the location and number of modal velocity peaks and the distribution of reflux. And has a modal velocity v _z ≈0.08mm s ^-1 Analog data set compared to MRI velocity profile of (c)The modal velocity peak of (2) is at a lower velocity v _z ≈0.02mm s ^-1 And has significantly more voxels associated therewith. Of course, MRI speed measurements are affected by experimental noise. Although it is estimated that in this case the noise-related overall speed error is relatively small, i.e. about 0.01mm s relative to the measured speed ^-1 But due to the lower SNR caused by the partial volume effect, larger noise-related errors may occur in the voxels at the grain edges. These errors in the velocity measurement can significantly reduce the sharpness of modal velocity peaks and generally lead to broadening of the velocity profile. With respect to reflux, the MRI dataset shows a greater degree of reflux than the reflux predicted by LBM extends further.

Voxel-by-voxel comparisons between the simulated velocity map and the experimental velocity map are shown in fig. 8. In general, good agreement between the two data sets is observed, as most voxels approach the straight reference line (white dashed line) with a slope of 1. For a relatively small number of voxels, a significant difference between MRI speed and analog speed is observed (note the logarithmic scale of the gray bars). FIG. 8 also shows that the number of voxels representing stagnant or nearly stagnant fluid in LBM simulations is shown to be significantly greater than v in MRI velocity maps _z ＝0mm s ^-1 Positive and negative speed ranges. This can be explained in part by experimental noise, especially at the pore-particle interface, where noise can impair the accuracy of the measured velocity due to partial volume effects.

The experimental and simulated velocity data are further analyzed using a pore scale method, wherein velocity information is obtained for each individual pore in the image, which is then compared between the two data sets. To achieve this, a (combined) binarized structural image is first generated from the separate binary masks of the MRI and LBM datasets by selecting only voxels that are shared between the two datasets. This makes it possible to later use the pore marking to identify the same pore consisting of the same number of voxels. Then use AVIZO ^TM The Chamfer distance map and marker-based watershed algorithm in (c) separated the pore space of the structural image into individual labeled pores. However, the method is thatThe labeled structural image and each velocity dataset (MRI or LBM) are then used individually as inputs in the analysis of the labeled pores to extract and correlate the flow properties of the individual pores.

The results of this analysis are shown in fig. 9A and 9B. FIG. 9A shows a pore-by-pore comparison between simulated velocity data and experimental velocity data; the dashed line in the curve represents a reference line with a slope of 1. While there is good overall agreement between the two data sets, the general trend is that LBM simulations slightly underestimate large pore speeds, but overestimate small and negative pore speeds; examination of fig. 6A-6C and 8 also reveals similar trends. However, this observation cannot be fully explained by experimental noise. A different view of the data is obtained in fig. 9B, where for both data sets, a fraction of the total flow is obtained Is plotted as a fraction of the number of those pores carrying such flow (N _p (i)/∑N _p ) Is a function of (2). As can be seen, the agreement between the experiment and the LBM simulation is excellent. An interesting aspect of the data shown in fig. 9B is to demonstrate that the flow field in Ketton rock exhibits significant non-uniformity in pore size, with 10% of the pores carrying ≡50% of the flow. These results are consistent with the visual inspection of the velocity field in fig. 6A-6C.

In the microscopic pores of Ketton, due to short T ₂ Relaxation time, it is difficult to capture information about the inhaled fluid by MRI. In water saturated Estaillades, due to longer T ₂ Relaxation times, this information can be more easily captured. To illustrate both the microporous and macroporous regions of fluid saturation in Estaillades rock, the magnetic resonance intensity image generated from complex-valued MRI data is co-registered with the μCT data of the rock; 2D (xy) slice images of co-registered MRI and μCT datasets are shown in FIG. 10, in which color-coded MR signal intensities represent water concentrations in microscopic and macropores, respectively. Fig. 10 shows that the microporosity occupies a large area of the rock, whereas only a few macroporosity can be identified in the rock structure, These macropores may carry a substantial portion of the flow through the rock. It can also be seen that there is no water present in the dense calcite particles of the Estaillades formation.

Fig. 11A-11I provide visual comparisons between LBM simulations of estailades and MRI velocity maps. Three different 2D slice images extracted from co-registered MRI and μct datasets are shown representing velocity components in the surface flow direction (z). Fig. 11A, 11D and 11G show 2D (xy) slice images of estailades obtained from 3D LBM simulation at a spatial resolution of 6 μm, while fig. 11B, 11E and 11H show co-registered images at a spatial resolution of 35 μm. Difference map (v) _z (MRI)-v _z (LBM)) is obtained by calculating the absolute difference between the coarse-grained LBM simulation and the MRI flow map. Three regions outlined with dashed boxes and labeled as regions a, b, c and d identify interesting flow patterns in the rock sample.

Considering the complex structure of the Estaillades formation, good qualitative consistency (neglecting the microporous region) was observed between the LBM and the three image slices of the MRI dataset, with some local deviations in some of the pores. Similar to Ketton, it can be seen in both LBM data and MRI data that the flow in estailades is highly non-uniform with several high-speed flow channels. MRI flowgrams also revealed that the fluid was stagnant or nearly stagnant in the microporous region of the rock (i.e., dark gray μct image intensity). Region "a" (shown in FIGS. 11A, 11B and 11C) represents macropores where LBM is predicted to have v _Z ≥5mm s ^-1 And the diameter extends to the formation of a high velocity flow path approximately half the pore width, while the MRI velocity profile shows lower velocities in the same pore. Interestingly, the acquired MRI velocity map indicates that two relatively small high-speed channels have been formed in region "a" rather than one, as predicted by LBM. In other macropores in this image slice with significant flow, LBM tends to underestimate flow rate, as is apparent from the difference plot (fig. 11C).

Region "b" (shown in fig. 11D, 11E and 11F) shows another high-velocity flow channel in the middle of the relatively large aperture. In this case, the LBM has accurately predicted the location and velocity of the flow field.

In the region "c" shown in fig. 11G, 11H and 11I, microporous rock particles can be seen. No flow information was obtained from this region using LBM simulation, as the pore size was below the imaging resolution of X-ray μct. However, with MRI, the average flow rate can be accurately measured at voxel sizes much larger than the pore size of the fluid saturated rock. In this case, the MRI flow chart shows that the microporous particles contain stagnant or nearly stagnant fluid.

Finally, the region "d" also shown in fig. 11G, 11H and 11I shows a complex flow behavior with adjacent positive and negative flow channels in a large macropore. This complex flow behavior has been predicted relatively accurately by LBM because of the good qualitative agreement observed in this region between MRI and LBM velocity data.

Visual inspection of fig. 11A-11I shows that while complex flow patterns are well predicted by LBM for some pores, significant differences between MRI data and LBM data can be seen for other pores.

To investigate these differences in more detail, a more quantitative analysis similar to that performed on Ketton was performed on a pore-by-pore basis. This analysis is performed on the macroporous region of the rock shared between the binarized structural images of the MRI dataset and the LBM dataset. The results of this analysis are shown in fig. 12A and 12B. This analysis produces a large number of small pores consisting of several voxels. Thus, for clarity, only pores with a radius greater than 35 μm are shown. Fig. 12A shows that there is a much greater dot spread and deviation from the reference line with a slope of 1 relative to the data observed for keton (fig. 9A). However, like the keton data, LBM simulations generally underestimate large pore speeds, but overestimate negative pore speeds. As can be seen in fig. 12A, there are a large number of positive and negative average pore z-velocities in the MRI velocity map with significant negative and positive average pore z-velocities Has pores of about 0mm s in LBM simulation ^-1 Is->Close examination of the LBM simulations shown in fig. 11A-11I also reveals similar trends, as the LBM simulations predict more highly localized flows than seen in MRI speed data. This observation is also supported by the curve shown in fig. 12B, where about 94% of the flow in the estailades rock is carried by only 10% of the pores based on LBM data. Similarly, data extracted from acquired MRI velocity maps indicate that flow through the macroporous region of the rock is highly heterogeneous, with 10% of the pores carrying ≡80% of the flow. Compared to Ketton (where 50% of the flow is carried by 10% of the pores), 60% more flow is passing through the same percentage of pores in the Estaillades formation, which is an indication of its highly heterogeneous rock microstructure.

Qualitative analysis of fig. 11A-11I shows that fluid saturated microporosity exhibits predominantly low flow rates. However, as is apparent from fig. 10, microporosity occupies a large portion of the total porosity in the estailades. In fact, microporosity may account for more than half of the total rock porosity in the Estaillades rock. To investigate the extent to which these fluid saturated micro-pores, which occupy a significant portion of the total porosity, affect the overall flow rate, a quantitative analysis was performed on MRI velocity maps whose velocities were encoded in the surface flow direction (Z). For this analysis, the flow rate Q for each axial (xy) image slice is calculated by summing the z-velocities of the individual voxels within a given slice and then multiplying by the cross-sectional area of the individual voxel _z . This analysis is performed on three separate images, namely a total porosity image, a macroporosity image, and a microporosity image. To generate Q for each image slice _z The velocity map is first multiplied by two masks that are generated by manual thresholding of the MRI intensity image (fig. 10).

A first mask is required to null background noise and separate the macroporous and microporous regions; in the context of this analysis, macroporous regions are defined as regions of large and medium size pores (i.e., pores that can be captured by μct approximately) that are colored green in fig. 10, but in the context of this study, microporous regions are defined as regions of much lower pore size than the image resolution—these regions are colored blue in fig. 10.

Partial volume effects need to be considered by a second mask calculated by averaging multiple binary masks at different thresholds, and in fact most voxels are only partially saturated with water in order to make the pair Q _z Is an accurate estimate of (a). Q (Q) _z The curves as a function of the position of the image slices of three different porosity images in the middle region of the rock are shown in fig. 13, which depicts the total porosity 72, macroporosity 74 and microporosity 76 of the rock. The dashed line represents the flow rate applied during the MRI experiment, i.e., at Q _z ＝0.03ml min ^-1 And (3) downwards.

It can be seen that the microporous region also contributes significantly to the total flow rate through the estailades plug, although most flows through the macropores. For example, for a rock region from about 0.5mm to 1.3mm along the z-direction, the macroporous and microporous regions have approximately the same contribution to the total flow through the rock. Relative to average total Q _z The average contribution of microporosity to flow in this region of the rock was estimated to be 36±5% (uncertainty estimated by changing the manual segmentation threshold by ±10% with respect to the selected threshold).

Note that in the Ketton MRI velocity image, the separate macropore network (as shown in fig. 6A-6C) accounts for the total Q in the rock _z Thus, microporosity is expected to contribute less to flow. This example shows that microporosity plays an important role in the fluid flow and transport process in microporous carbonates.

The method may be performed in any order and repeated as necessary. Part or all of the method may be performed.

While embodiments have been described with reference to various embodiments and modes of use, it should be understood that these embodiments are illustrative and that the scope of the inventive subject matter is not limited in this respect. Many variations, modifications, additions, and improvements are possible. For example, one or more images may be performed using one or more of the techniques herein. Various combinations of the techniques provided herein may be used.

Multiple examples may be provided for components, operations, or structures described herein as a single example. Generally, structures and functionality presented as separate components in the exemplary configurations may be implemented as a combined structure or component. Similarly, structures and functionality presented as a single component may be implemented as separate components. These and other variations, modifications, additions, and improvements may fall within the scope of the inventive subject matter.

Claims (16)

1. A method for calibrating direct flow simulation of a rock sample, the method comprising the steps of:

-providing a 3D image of a rock sample, wherein the 3D image is composed of a plurality of voxels;

-generating a segmented structure image of the rock sample from the 3D image by selecting each voxel of the plurality of voxels to represent a void space or a solid material;

-simulating fluid flow on the segmented structure image with a direct flow simulation;

-generating a 3D spatially resolved fluid velocity map of one or more fluid phases at a pore scale resolution using pulsed field gradient magnetic resonance imaging; and

-comparing the simulated fluid flow and the 3D spatially resolved fluid velocity map;

thereby calibrating the direct flow simulation on the rock sample.

2. The method of claim 1, wherein the step of comparing comprises co-registering the simulated fluid flow and the 3D spatially resolved fluid velocity map.

3. The method of claim 1, wherein the direct flow simulation is based on the lattice boltzmann method.

4. The method of claim 1, wherein the pulsed field gradient nuclear magnetic resonance is combined with a RARE MRI pulse sequence.

5. The method of claim 1, wherein the step of generating a 3D spatially resolved fluid velocity map comprises the steps of: spatially resolving a fluid flow in the rock sample; reconstructing a spatially resolved fluid flow encoded phase difference map for each of the one or more fluid phases in a flow state and a static state using compressive sensing; and calculating a flow velocity map for each of the fluid phases.

6. The method of claim 5, wherein the step of spatially resolving the fluid flow is performed by velocity mapping to provide an average velocity of the fluid flow for each voxel within the 3D image.

7. The method of claim 6, wherein the velocity mapping is performed with a combination of pulsed gradient spin echo sequences and MRI pulse sequences.

8. The method of claim 7, wherein the MRI pulse sequence is a rareiri pulse sequence.

9. The method of claim 5, wherein the step of spatially resolving the fluid flow is performed using spatially resolved displacement propagators to provide a fluid flow distribution for each voxel within the 3D image for a predetermined time.

10. The method of claim 9, wherein the spatially resolved shift propagators are acquired by a combination of alternating pulse gradient stimulated echo sequences and MRI pulse sequences.

11. The method of claim 10, wherein the MRI sequence is a rareiri pulse sequence.

12. The method of claim 1, wherein the 3D image of the rock is obtained by X-ray computed tomography.

13. The method of claim 1, wherein the step of comparing the simulated fluid flow and the 3D spatially resolved fluid velocity map comprises quantifying a contribution of microporosity in the rock sample to fluid flow.

14. The method of claim 1, wherein the calibrated direct flow simulation is used to make decisions regarding recovery of hydrocarbons from a subsurface formation.

15. The method of claim 1, wherein the calibrated direct flow simulation is used to make decisions regarding carbon capture and sequestration in the subsurface formation.

16. The method of claim 1, wherein the calibrated direct flow simulation is used to make decisions regarding geothermal extraction processes.