WO2019151889A1

WO2019151889A1 - A method for determining a three-dimensional spatial distribution of porosity in a sample of a heterogeneous porous medium

Info

Publication number: WO2019151889A1
Application number: PCT/RU2018/000058
Authority: WO
Inventors: Igor Andreevich VARFOLOMEEV; Dmitry Alexandrovich Korobkov; Ivan Victorovich YAKIMCHUK
Original assignee: Schlumberger Technology Corporation; Schlumberger Canada Limited; Services Petroliers Schlumberger; Schlumberger Technology B.V.
Priority date: 2018-02-02
Filing date: 2018-02-02
Publication date: 2019-08-08
Also published as: RU2020128473A3; RU2020128473A

Abstract

A method for obtaining 3D spatial distribution of porosity (porosity map or porous model) inside a sample of a heterogeneous porous medium comprises obtaining a 3D microstructural image of the sample by a 3D microstructural imaging procedure and measuring a total porosity of the sample. Then, values in the obtained 3D microstructural image that correspond to fully void and fully solid voxels are determined and the obtained 3D microstructural image data is normalized using the determined values corresponding to the fully void and the fully solid voxels. Finally, a 3D spatial distribution of porosity of the sample is created by a computing device using the normalized 3D microstructural image and the measured total porosity.

Description

A METHOD FOR DETERMINING A THREE-DIMENSIONAL SPATIAL DISTRIBUTION OF POROSITY IN A SAMPLE OF A HETEROGENEOUS POROUS MEDIUM

Field of the disclosure

The disclosure relates generally to the field of estimating properties of heterogeneous porous media. More specifically, the invention relates to methods for imaging and porosity analysis of a sample of a heterogeneous porous medium.

Background

Porous medium is an important type of matter, which can be found almost everywhere around us. Reservoir rocks in oil and gas collects and filtrates hydrocarbons through their pores. Properly designed porous structure in fabric makes modem clothes more effective when interacting with moisture. Porosity and exact 3D geometry of pore space is an extremely critical factor in performance of construction materials (building foundations, bridge supports, etc.). These and lots of other examples vividly demonstrate the importance of ability to measure and analyze porosity of a sample and even three-dimensional spatial distribution of porosity within a sample.

Depending on a sample size and scale of pores different measuring systems exist. For example, in oil and gas industry, the logging is widely used for getting information about formation porosity. However, well logs cannot be directly interpreted without “calibration” to lab measurements. In laboratory, core plug porosity can be directly measured with a special device (via gas saturation) at ambient conditions as well as at high pressure and temperature. However, these and most of other methods evaluate only single integral value, e.g. total porosity.

Nowadays, image-based approach for studying sample properties become a well- recognized alternative or supplement to conventional laboratory techniques. Thus, 3D digital models of porous samples (e.g. rocks) with consequent numerical simulation of physical properties are important all-sufficient results, which start playing crucial role in modem industrial workflows

(see, for example, E. Maire and P. J. Withers. Quantitative X-ray Tomography. International Materials Reviews. 2014; 59(1 ): 1-43).

The success of application of such technologies including mentioned above depends significantly on the quality of input data - digital models. Existing solutions for creating 3D digital models are based on 3D imaging of the sample. The geometrical accuracy of achieved 3D digital models is limited by spatial resolution of 3D imaging technique. Commonly X-ray microtomography is used as a 3D imaging method with best spatial resolution - 1 pm. At the same time, it is well known that pore size distribution can spread from nanometers to millimeters depending on a sample nature. As a results, pore space is partially lost (not resolved) in 3D digital models (e.g. digital rock models).

Summary

Disclosed is a method for obtaining 3D spatial distribution of porosity (porosity map or porous model) inside a sample of a heterogeneous porous medium. The method comprises obtaining a 3D microstructural image of the sample by a 3D microstructural imaging procedure and measuring a total porosity of the sample. Then, values in the obtained 3D microstructural image that correspond to fully void and fully solid voxels are determined and the obtained 3D microstructural image data is normalized using the determined values corresponding to the fully void and the fully solid voxels. Finally, a 3D spatial distribution of porosity of the sample is created by a computing device using the normalized 3D microstructural image and the measured total porosity.

Brief description of the drawings

Features and advantages of the described implementations can be more readily understood by reference to the following description taken in conjunction with the accompanying drawings.

Fig. 1 is a flowchart illustrating an example method in accordance with some embodiments;

Fig. 2 shows an example of 2D slice of an obtained 3D microstructural image (x-ray microtomography) of a porous sample (carbonate rock);

Fig. 3A-3F illustrates the problem of limited spatial resolution of 3D microstrucural imaging: Fig. 3A shows an artificial structure without partially void/solid voxels (a monomineral sample); Fig.3B - a synthetic microstructural image of the artificial structure of Fig.3 A with insufficient spatial resolution; Fig. 3C - the Fig,3B image with noise; Fig. 3D - the histogram of Fig. 3A; Fig.3D - the histogram of Fig. 3B; Fig. 3F - the histogram of Fig. 3C;

Fig. 4A-4F demonstrates a special preprocessing step for polymineral samples: Fig. 4A demonstrates an artificial structure without partially void/solid voxels and with two mineral types; Fig.4B shows a synthetic noisy microstructural image with resolution limitations; Fig.4C shows the result of processing of the image shown on Fig. 4B; Fig.4D is the histogram of Fig. 4A; Fig.4E - the histogram of Fig.4B; Fig. 4F - the histogram of Fig.4C;

Fig. 5A-5B helps with understanding of final creation of 3D spatial distribution of porosity (3D porosity map); Fig. 5 A is the histogram of the image that corresponds to the condition of linear dependence between voxel values and sample local density in the voxel volume, values P_\ and Po are average voxel values for fully void and fully solid voxels correspondingly; Fig. 5B shows how the final 3D spatial distribution of porosity is obtained through determination of optimal Pi or Po for satisfying the condition of equal total <t>r and digital Fn porosities of the sample.

Fig. 6 illustrates a computing system in accordance with some embodiments;

Fig. 7A-7E shows an example case of constructing 3D porous model of carbonate sample and further calculation of its electrical formation factor F: Fig. 7A shows a part of 2D slice of an initial 3D microstructural image (X-ray microtomography) of the carbonate sample; Fig. 7B| shows the result of automated thresholding procedure, Fig. 7C shows the result of manual thresholding, which satisfies the condition of equal total F_G and digital Få porosities of the whole sample; Fig. 7D shows the obtained with disclosed invention 3D spatial distribution of porosity; Fig. 7E shows the porosity F and formation factor F values numerically calculated for achieved 3D spatial distributions of porosity Fig. 7B-7D.

Fig. 8A-8C illustrates the preservation of connectivity between pores after application of the disclosed method in comparison with conventional thresholding approach: Fig. 8A shows a

small 2D region of original 3D microstructural image; Fig. 8B shows a result of thresholding and Fig. 8C shows a corresponding part of created 3D spatial distribution of porosity.

Detailed Description

The disclosure is aimed on significant decrease of a negative impact of not resolved porous space, thus reducing a probability of making mistakes in further business decisions based on 3D porosity distribution of a sample. The invention expands the scope of applicability of 3D imaging techniques, such as but not limiting to 3D X-ray (micro)tomography. Moreover, in many applications increase of imaging resolution can be too expensive and, thus, less profitable than utilization of disclosed invention.

Fig. 1 shows a flowchart in accordance with one or more embodiments.

Block 1 corresponds to a heterogeneous media sample. For example, it can be a core sample consisting of rock minerals that was extracted from a near wellbore area or a concrete sample.

In Block 2 an initial 3D microstructural image (shown in Block 3) of the sample is obtained via 3D microstructural imaging procedure. The 3D microstructural image shown in Block 3 is an image that explicitly or implicitly reflects internal geometry and mineral heterogeneity in the sample with limitations related to spatial and mineral resolutions.

Digital representation of a 3D image is a 3D array of scalar or vector values. Each element of that array corresponds to a point in the image (pixel for 2D image and voxel for 3D image). A large number of methods and corresponding devices can provide a 3D microstructural image. Among them are such methods as X-ray micro- (nano-) tomography, X-ray fluorescence microtomography, neutron microtomography, 3D FIB-SEM, etc. and well-known devices for their realization.

It is important that applicable 3D imaging methods should satisfy the following requirement: values in voxels should be in a reversible function dependence on substance density (therefore, porosity) for a homogeneous substance. In such a case, the values in voxels can be recalculated to linearly dependent on homogeneous substance density. X-ray microtomography is a good example of technique satisfying this requirement (NIST: X-Ray Mass Attenuation

Coefficients - Section 2 - http://phvsics.nist.gov/PhvsRefData/XravMassCoefichap2.htmn. An example of X-ray microtomography image of a carbonate is presented on Fig. 2, where greyscale of the image from black to white corresponds to sample local X-ray attenuation m from low to high absorption that linearly depends on density for particular substance.

It should be noted that in general case every imaging device has an apparatus function of relation between value in the point of image and corresponding physical value in this point. Typically, this function is more or less known or can be estimated and, thus, the image can be converted to at least linear correspondence between the value in the point of an image and corresponding physical value in this point.

For convenience, hereinafter a 3D microstructural image is an image already recalculated to linear dependence on density in the voxel so, that increase of density leads to increase of value in the image.

The next step is shown in Block 4 in Fig. 1 , where the heterogeneous media sample is studied in a laboratory by any of numerous well-proved porosity measurement techniques. In one or more embodiments, it can be gas porosimetry. It is believed, that such techniques determine total porosity <Dr (shown in Block 5) of the sample in contrast to 3D microstructural image processing, which can partially lose the volume of void space due to spatial resolution limitations.

In case of dealing with core plugs the porosity value can be measured by direct methods - by gas using Boyle’s Law Single Cell Method or by saturation of the sample with a liquid (by brine, oil or kerosene) i.e. filling the pore space with a fluid and measurement of the fluid amount by weighting. In another embodiment of invention, effective porosity can be measured by mercury porosimetry. In some cases, information on porosity can be obtained from log data.

In case of multicomponent (multiphase) sample (e.g. polymineral rock, multilayer structure) an additional 3D mapping of chemically homogeneous components is recommended. As an example, for polymineral rock samples 3D mineral mapping (Block 6) is required. Hereinafter a mineral is a specific chemical substance (carbon, gold, specific metal alloy, etc ), or a specific mineral type (quartz, calcite, pyrite, etc.), or a specific mixture of minerals (e.g. 30% of feldspar and 70% of quartz), or a specific mixture of chemical substances (e.g. different layers in a multilayer shell of a sample). Any of the above is called as mineral all over the invention

description. Samples with multimineral (with just determined notation of“mineral”) composition are referred as polymineral in the whole description. Hereinafter a 3D mineral map is an image where every point of the image can be interpreted as a particular mineral. As a result, a 3D mineral map (Block 7) of the sample is obtained, which corresponds to 3D distribution of chemically homogeneous components, in general case.

In one or more embodiments, 3D mineral mapping (Block 6) can be a synchrotron microtomography with monochromatic X-ray beam (F. Fusseis, X. Xiao, C. Schrank and F. De Carlo, "A brief guide to synchrotron radiation-based microtomography in (structural) geology and rock mechanics," Journal of Structural Geology, vol. 65, pp. 1-16, 2014). The ability to resolve different minerals can be increased by scanning the sample several times with various energies (wavelengths) of X-ray beam - multi-energy microtomography. In another embodiment, SEM complemented with Focused Ion Beam (FIB) system can be used. It allows etching and imaging the surface of the object under study slice by slice. This technique may provide 3D image of the subsurface volume with ~l 0 nm resolution and ~l 0 pm field of view. Backscattered electron (BSE) detection together with EDS can provide mineral information for acquired 3D image. In one or more embodiments, other 3D imaging techniques sensitive to chemical and mineral content can be applied. E.g., 3D X-ray fluorescence microtomography, 3D confocal Raman spectroscopy, 3D X- ray topo-tomography.

In another embodiments, modem indirect methods can be applied. E.g. gray values in 3D X-ray microtomography (microCT) image can be linked with different densities and chemical composition, such as different mineral (WO2013058672). The process can be improved by correlating gray values of 3D microCT image with spatially registered 2D mineral distribution obtained with scanning electron microscopy and energy dispersive spectroscopy (US20150104078A1). More sophisticated approach consists in extracting (calculating) various local features in points of 3D microCT image prior to matching with 2D mineral maps image (I. Varfolomeev, I. Yakimchuk and B. Sharchilev, "Segmentation of 3D image of a rock sample supervised by 2D mineralogical image," // 2015 3rd IAPR Asian Conference on Pattern Recognition (ACPR), Kuala Lumpur, 2015, pp. 346-350. doi: 10.1 109/ACPR.

In Block 8 a special transformation of 3D microstructural image is performed - 3D porosity mapping. As a result, a 3D spatial distribution of porosity or 3D porosity map (Block 9) is constructed, i.e. a 3D image f ( c,g, z) containing in each voxel the fraction of pore space in the voxel volume (intrinsic porosity of voxel). The value of“1” would be assigned to a fully void voxel, while a totally solid voxel would have the value of“0”. In contrast to classical binarization of 3D microstructural image, voxels with intermediate values between“0” and“1” also exist in porosity map/distribution f(c,g, z) . Such values f represent intrinsic porosity of partially void/solid voxels.

One important note should be made here. A totally void voxel means that there is no mineral (solid) phase in the corresponding volume element of the sample, while it is not necessarily physically completely empty (vacuum). Normally, it is filled with air, water, oil or any other types of fluids or gases. Hereinafter, such voxels would be called as totally void voxels with intrinsic porosity f = 1.

The 3D spatial distribution of porosity is achieved by accounting during segmentation the laboratory measured total porosity (Block 5) and additional accounting the 3D distribution of minerals (chemically homogeneous components) in case of polymineral (multicomponent) sample. The general idea of the procedure is based on the following assumptions:

a) Most of the voxels of the 3D microstructural image are intrinsically monomineral, i.e. each voxel contains only voids and some single type of a mineral (material, chemically homogeneous substance) inside, not a mixture of minerals. The type of mineral does not necessarily the same for all voxels.

b) According to above-mentioned requirement to imaging method and (a), values in such voxels of the 3D microstructural image almost linearly decrease with increase of intrinsic porosity of voxel for every type of mineral (material, chemically homogeneous substance) independently.

Thus, in case of a monomineral sample an ideal (no noise or any other artefacts) 3D microstructural image (artificial example demonstrated in Fig. 3B) should contain values in voxels in the range from I_\ to Io (Fig. 3E - image histogram of Fig. 3B), where h corresponds to totally void voxels (intrinsic porosity f = 1) and Io - to fully solid voxels (intrinsic porosity f = 0). True artificial structure (without partially void/solid voxels) that corresponds to this image is shown in

Fig. 3A, where dark grey round pores are located in light grey mineral matrix. It has only two values h and 7o (Fig. 3D - image histogram of Fig. 3 A). Real images are always prone to noise and other artifacts (Fig. 3C). As a result, values lower than 7i and higher than 7o can be found in the 3D microstructural image (Fig. 3F - image histogram of Fig. 3C). In one or more embodiments, values 7o (and 7i) can be determined manually by averaging the group of voxels belonging to internal area of resolved totally void pore (or totally solid mineral grain). In another embodiment, the average value can be replaced by statistical mode (the most frequent value).

Consequently, in case of a monomineral sample in accordance to above stated assumptions, image values in the range from 7i to 7o can be easily converted to intrinsic porosities.

In case of a polymineral sample, previous concept is directly valid only for a set of monomineral voxels of each mineral type M independently. More precisely, the value 7i in all totally void voxels is the same (within noise level) for every mineral type M, but the 7o( ) is the function of mineral type M. In other words, if the sample is a dolomite with anhydrite inclusions (only two minerals are considered for simplification) all voxels can be classified as‘porous dolomite’ or‘porous anhydrite. The values 7o(A7) still can be defined by the same way as in monomineral case. This situation is illustrated in Fig. 4. Fig. 4A demonstrates artificial structure without partially void/solid voxels and with two mineral types - anhydrite at the top half of the image and dolomite at the bottom. Fig. 4D is the histogram of Fig. 4A. As one can see, three values are observed in the image: h (corresponds to void voxels), 7OA (corresponds to dolomite voxels) and 7OB (corresponds to anhydrite voxels). Real noisy image with resolution limitations is shown in Fig.4B. Corresponding histogram is presented in Fig. 4E.

Thus, a special processing of such 3D image is required for making the values in the image linearly dependent on intrinsic densities. In one or more embodiments, for transforming polymineral case to monomineral-like, one can process each non-totally void voxel of the image in the following manner:

where P(x,y,z ) - is a value in the voxel

after processing of the 3D microstructural image /(x,_y,z), M(x,y,z) - is the 3D mineral map already derived by 3D Mineral Mapping procedure, m - is the estimation of the physical value imaged by selected technique in fully void voxels, and U_{M(x y} - the same for fully solid voxels that corresponds to mineral type M. The method for obtaining these estimations strictly depends on the method of 3D microstructural imaging.

In one or more embodiments, in case of 3D X-ray microtomography such values can be calculated theoretically (NIST: X-Ray Mass Attenuation Coefficients - Section 2 - http://phvsics.nist.gov/PhvsRefData/XravMassCoef/chap2.htmft. In another embodiment, 3D microstructural imaging of mineral (gas, fluid) standards would be a universal solution for every imaging method. In other embodiments of the invention, determination of statistical mode value in some neighborhood of every voxel (x,y,z) can be treated as P_{M[x y :} . The size of neighborhood should be large enough to provide predominance of totally solid voxels in it.

Having processed 3D microstructural image, one will achieve a 3D image P(x,y, z ) with values linearly depending on intrinsic porosities - smaller value P(x, y, z ) corresponds to a larger value of intrinsic porosity of the voxel (x, y, z) and vice versa. Similarly to /o and /i, let Po and P be the values of P{x,y, z) that represent fully solid and fully void voxels correspondingly.

According to (1 ), P is equal to 0, while P₀

. For correct 3D microstructural imaging

— M pore

and proper estimation of m_{rq ί} and m_{M[c y z)} the value of Po is the same (within noise level) for every mineral type.

Returning to artificially created polymineral image presented in Fig. 4, one can see in Fig. 4C the result of transformation (1) of the Fig. 4B. After processing, the image contains values in the range from P to Po (Fig. 4F - image histogram of Fig. 4C) that linearly depend on intrinsic porosities, as in case of monomineral sample.

For convenience, both cases (monomineral and polymineral) can be further presented in one notation. Let P{x,y,z ) = /(« :) in a monomineral case, then Po and P would be just equal to 7o and h. Thereby, all considerations and statements below would be correct for both mono- and polymineral cases.

Finalizing and concluding the description of notions To and 7i, h(M), P{x,y,z), Po and P , the next statement can be claimed. To proceed with construction of 3D spatial distribution of porosity after obtaining 3D microstructural image and measuring total porosity of the sample one need to determine values in the obtained 3D microstructural image that correspond to fully void and fully solid voxels for every mineral type in the sample, i.e. 7i and 7o (or 7o(A7) in case of a polymineral sample).

In general (mono- and polymineral) case, the next required step is normalization of the 3D microstructural image data using the determined values that correspond to the fully void and the fully solid voxels, i.e. h and 7o (or 7o(A7) in case of a polymineral sample). The normalization procedure consists in calculating P{x,y,z) . As mentioned above, for the monomineral case

P(x,y,z) = I (x ,y,z ) , while the polymineral case is described by equation (1).

Finally, one can create, by a computing device, a 3D spatial distribution of porosity of the sample using the normalized 3D microstructural image P{x,y, z) and the measured total porosity. This procedure is described in detail below.

All voxels smaller than i are considered as fully void. They relates to spatially resolved porous volume VRP. As it is mentioned above, usually this volume VRP is smaller than actual total porous volume of the sample VT. Difference between these values VNP = (VT - VRP) is allocated to a set of (micro )porous regions non-resolved by 3D microstructural imaging. All voxels larger than Po are considered as fully solid. They relate to spatially resolved solid mineral matrix volume VRS. Similarly, VNS - is volume of a solid mineral matrix non-resolved by 3D microstructural imaging. The sum of VRS, VRP, VNS and VNP is equal to a volume of the imaged sample V_å . For convenience, hereinafter all volume values have voxels as units. Thus, VRS is a number of totally solid voxels (f = 0), V_RP is a number of totally void voxels (cp = 1), and V_N = V_NS + VNP corresponds to a number of voxels that are partially void/solid.

Finally, to derive 3D distribution of intrinsic voxel porosities p{x,y,z) the correspondence between values P{x,y,z) and porosities should be defined. In one or more embodiments, it can be done in accordance to previously determined values Po and P_\ by using the following formula:

Dash line in Fig. 5A illustrates the behavior of cp{P(x,y,z) . This line is overlaid on histogram of image Fig. 4C.

After processing according to (2), the whole 3D microstructural image, the overall image

porosity 0_å = ^{ix,y z)} - can differ from experimentally measured total porosity Or . In cases p_å

where there is a confidence that both values should be equal (close within some range) to each other, another embodiment of the disclosure is proposed - the value of Po can be found from the condition of equality between the experimental total porosity Or and the derived from image overall porosity O_{S .} This approach can be reasonable if it is hard to determine o just from the image, e.g. all voxels are partially porous. The procedure can be organized in different ways. In one or more embodiments, straightforward iterative approach can be used (Fig. 5B). Changing incrementally the value of Po from Pi to larger values, will monotonically increase the 0_å from

— to 1 (for P₀ ® +ao ). For correctly performed total porosity measurement and 3D

På

V V

microstructural imaging total porosity Or satisfies the following condition:—— < O_r < 1—— .

p_å på

Thus, one can iteratively reach the satisfaction of condition F_S = Or with a reasonable accuracy. In other embodiments of the invention, alternative more sophisticated approaches can also be

applied (https://en. wikipedia. org/w/index.php?title=Root-fmding_algorithm&oldid=696250l 66) for determining optimal Po and, consequently, 3D distribution of intrinsic voxel porosities

H^x>y>^z) ·

It should be noted, that in case of absence of totally void and/or solid voxels in the obtained 3D microstructural image the proper values of 7i and/or 7o, i and/or Po will be outside the range of values observed in the image, while VRP = 0 and/or VRS = 0, correspondingly. Thus, the value of 7i cannot be determined from the image and the procedure (1) cannot be performed. Although such situation is rather unlikely and extraordinary, still it can be treated by the disclosre if at least one of the values 7i, (7_o(M)} can be obtained from the image.

In such specific case, for implementation of the disclosure in a monomineral case (procedure (1) is not performed), described above solution works without any changes if h (and consequently Pi) can be obtained from the image (totally void voxels can be found in the image). Then Po can be determined as stated above. If only 7_o (and consequently Po) can be obtained from the image (totally solid voxels can be found in the image), the process should be changed slightly. Similarly to above discussed embodiments, instead of Po the value of Pi can be found from the condition of equality between the experimental total porosity F₇ and the derived from the image overall porosity F_{S .} The procedure can be organized in different ways as well. In one or more embodiments, straightforward iterative approach can be used. Decreasing the value of Pi from Po

V V V

to lower values, F_S will monotonically decrease from 1—— to 0. As far as— ¹— < F_G < 1—— ,

one can iteratively reach the satisfaction of condition F_S = F_G with a reasonable accuracy. In other embodiments of the invention, alternative more sophisticated approaches can also be applied for determining optimal Pi and, consequently, 3D distribution of intrinsic voxel porosities f(c, y,z) .

For implementation of the disclosure in a polymineral case, if 7i (and consequently Pi) can be obtained from the image (totally void voxels can be found in the image) initially described solution also works without any changes. Procedures (1) and (2) with determination ofPo from the condition of equality between the experimental total porosity F_G and the derived from the image overall porosity F_S are absolutely the same. In other embodiments, when 7i (and consequently Pi)

cannot be defined directly from the image (all voxels are partially solid), but at least one of {7o( )} can be obtained - I_ref, procedure (1) and (2) should be replaced with more general optimization problem:

Here, (3.1) reflects the fact of linear dependence of image values I{x,y,z) on effective physical values m_naa1 (x,y, z) imaged by selected technique, (3.2) is the relation (3.1) applied to the totally solid voxel of a reference mineral, image value of which can be obtained. Expression (3.3) defines the effective physical value of m_n¥c1 (x,y,z) for known physical value of pore space m_rok , solid space m_{M{c y ,,} and intrinsic porosity <p{x, y, å) in this voxel (x,y,z) . (3.4) is the expression for 3D distribution of intrinsic porosities with unknown slope coefficient A. The optimization problem (3.5) consists in determining proper value of A that provides equality of

experimental total porosity F_G and overall image porosity F_S = -

Actually, (3.5) is just an equation for A with a big number of terms (for each voxel) similar to (3.4). It can be easily solved by any of well-known methods, including already mentioned iterative approach.

In one or more embodiments, the issue of absence of totally void and/or solid voxels in the object can be fixed by artificial addition of totally void and/or solid object in it and the whole workflow of the invention should be repeated. Besides, such approach is the only method in case of impossibility to derive any of the values 7i, (7o(M)} from the image.

Obtained 3D distribution of intrinsic voxel porosities f (c,g,z) is a 3D digital model of the sample with mapped porosity in each voxel or 3D porosity map, i.e. desired 3D spatial distribution of porosity in the sample.

In one or more embodiment, the construction of final 3D spatial distribution of porosity in the sample can include additional post-processing. For example, in real practice, direct 3D porosity mapping (as described above) may lead to formation of partially solid voxels (<p < 1 ) surrounded by fully void voxels (f = 1), which contradicts to the real physics. Alternatively, partially void voxels can be surrounded by fully solid voxels and, thus, do not have any impact on a process under study, e.g. fluid or electrical current flow. In one or more embodiments, such isolated partially solid voxels in fully void space (or vice versa) can be replaced by the surrounding phase (void or solid).

In one or more embodiments, all necessary calculations and image processing required for constructing 3D spatial distribution of porosity in the sample are performed in a parallel mode.

As a result, 3D porosity map (3D spatial distribution of porosity in the sample) is constructed (Block 9 on Fig.1 ). In one or more embodiments, it can be stored as a 3D array of intrinsic porosity values in each voxel. This 3D dataset can be stored as one file or a set of 2D slices.

The obtained 3D spatial distribution of porosity in the sample characterizes the sample in terms of structure peculiarities. Besides, such distribution can be used as a digital model of the heterogeneous media for numerical simulations of various physical phenomena: gas/fluid flow, electromagnetic, mechanical strength, etc.

A system for creating a 3D spatial distribution of porosity in the sample of a heterogeneous media comprises an image-producing device configured to produce an initial 3D microstructural image of a sample. This device is selected from a group of devices providing such methods as X- ray micro- (nano-) tomography, X-ray fluorescence microtomography, neutron microtomography, 3D FIB-SEM, etc.

The method requires using a computing device coupled to the image-producing device. The computing device may include hardware, software, firmware, or a combination thereof. Various components of the computing device are described below with reference to Fig. 6.

As shown in Fig.6, the computing device may be of virtually any type regardless of the platform being used. For example, the computing device may be one or more mobile devices (e.g., laptop computer, smartphone, smartwatch, personal digital assistant, tablet computer, or other mobile device), desktop computers, servers, blades in a server chassis, or any other type of computing device or devices that includes at least the minimum processing power, memory, and input and output device(s) to perform one or more embodiments of the invention. For example, as shown in Fig. 6, the computing device may include one or more computer processor(s) 14, associated memory 15 (e.g., random access memory (RAM), cache memory, flash memory, etc.), one or more storage device(s) 16 (e.g., a hard disk, an optical drive such as a compact disk (CD) drive or digital versatile disk (DVD) drive, a flash memory stick, etc.), and numerous other elements and functionalities. The computer processors) 14 may be an integrated circuit for processing instructions. For example, the computer processor(s) may be one or more cores, or micro-cores of a processor. The computing device may also include one or more input device(s) 17, such as a touchscreen, keyboard, mouse, microphone, touchpad, electronic pen, or any other type of input device. Further, the computing device may include one or more output device(s) 18, such as a screen (e.g., a liquid crystal display (LCD), a plasma display, cathode ray tube (CRT) monitor, e-ink display, projector, or other display device), a printer, external storage, or any other output device. One or more of the output device(s) may be the same or different from the input device(s). The computing device may be connected to a network 19 (e.g., a local area network (LAN), a wide area network (WAN) such as the Internet, mobile network, or any other type of network) via a network interface connection. The input 17 and output 18 device(s) may be locally or remotely (e.g., via the network 19) connected to the computer processors) 14, memory 15, and storage device(s) 16. Many different types of computing devices exist, and the aforementioned input and output device(s) may take other forms.

Software instructions in the form of computer readable program code to perform one or more embodiments may be stored, in whole or in part, temporarily or permanently, on a non- transitory computer readable medium such as a CD, DVD, storage device, a diskette, a tape, flash memory, physical memory, or any other computer readable storage medium. Specifically, the

software instructions may correspond to computer readable program code that when executed by a processor(s), is configured to perform one or more embodiments of the method.

Further, one or more elements of the aforementioned computing device may be located at a remote location and connected to the other elements over a network (19). Further, embodiments may be implemented on a distributed system having multiple nodes, where each portion of an embodiment may be located on a different node within the distributed system. In one or more embodiments, the node corresponds to a distinct computing device. Alternatively, the node may correspond to a computer processor with associated physical memory or to a computer processor or micro-core of a computer processor with shared memory and/or resources.

Example of invention implementation

A part of 2D slice of initial 3D microstructural image (X-ray microtomography) of a carbonate sample is presented in Fig. 7A. Laboratory measured total porosity Ft is equal to 17.2%. Standard thresholding technique provides only 9.2% of porosity (Fig. 7B). To achieve experimentally measured porosity one should select excessively high threshold value (Fig. 7C).

By applying disclosed invention with defined from the image Pi and numerically optimized value of o, 3D spatial distribution of porosity in the sample was constructed (Fig. 7D).

To emphasize the importance of the disclosed method and its result an electrical field distribution and corresponding electrical conductivity (or formation factor F) for this rock sample (saturated with a particular brine) was evaluated numerically.

This simulation requires the following:

(a). Conversion of constructed 3D spatial distribution of porosity in the sample to 3D distribution of local electrical conductivity a(x,y,z ) - conductivity of every voxel.

(b). Determine electrical field distribution under given boundary conditions.

(c). Calculate effective conductivity characteristics of the whole digital model.

Step (a) can be performed using the following rule:

- for totally solid voxels conductivity s ( , y, z) = cr_{M[x y z)} ;

- for totally void voxels conductivity s ( x, y, z) = cr_hnnc ;

for partially solid/void voxels one possible way to assign conductivity is

Here M{x, y, z ) is a mineral type in voxel with coordinates (x,y, z) , s_{M{c y} _ is the corresponding electrical conductivity, a_bnne is the conductivity of brine. In this particular example, the sample was a dolomite with anhydrite inclusions and the following values have been selected:

°bri«e = 2.5 S/m, a_dolomle = 0.3 x l0 ⁴ S/m, a_®dnte = 10 ⁹ S/m.

Numerical simulation of electrical field distribution can be performed by iteratively solving Laplace equation for electrical potential value using obtained conductivity distribution with specific boundary conditions (Patankar S.V,“Numerical Heat transfer and fluid flow”, 1980, pp.59-66). Boundary conditions can be set in such a way to correspond the considered physical problem.

For example, for a cubic model the following boundary conditions can be used:

constant electric potential difference applied to sides perpendicular to the direction of electric field;

- zero normal derivative of the potential on the other boundaries.

Other numerical approaches are also possible, e.g., based on inversion ofKirchhoffs circuit laws matrix.

Using electric potential and conductivity distributions in media electrical currents distribution can be obtained (in accordance with the Ohm’s law). Effective conductivity, formation factor and other conductivity characteristics can be determined by means of averaging conductivity currents distribution for a given applied electrical potential difference (https://en.wikipedia.org/w/index.php?title=Archie%27s_law&oldid=692860l88).

As a result (Fig. 7E), for traditional thresholding approach (column B) numerical formation factor F was -100 times higher (Ft, = 3051) than should be (Exp_eriment = 30) for such rock (known form independent experimental studies). For digital model created with overestimated threshold (column C), formation factor is still ~5 times higher (F_c = 160). In contrast, disclosed invention (column D) provides good agreement between numerical (E_d = 26) and experimental formation factors (X_experiment = 30). One of the main reasons for this consists in preserved connectivity of pore

space in the sample. Fig. 8 illustrates this feature: small 2D region of original 3D microstructural image is shown in Fig. 8A, result of thresholding - Fig. 8B and corresponding part of 3D spatial distribution of porosity in the sample - Fig. 8C.

Claims

1. A method of determining a three-dimensional spatial distribution of porosity in a sample of a heterogeneous medium, the method comprising:

- obtaining a 3D microstructural image of the sample by a 3D microstructural imaging procedure,

- measuring a total porosity of the sample,

determining values in the obtained 3D microstructural image that correspond to fully void and/or fully solid voxels,

normalization of the obtained 3D microstructural image data using die determined values corresponding to the fully void and the fully solid voxels,

creating, by a computing device, a 3D spatial distribution of porosity in the sample using the normalized 3D microstructural image and the measured total porosity.

2. The method of claim 1 wherein the sample of the heterogeneous medium is a core.

3. The method of claim 1 wherein the 3D microstructural image is obtained by X-ray micro-

(nano-) computed tomography.

4. The method of claim 1 wherein the 3D microstructural image is obtained by neutron microtomography.

5. The method of claim 1 wherein the obtained 3D microstructural image is a vector image.

6. The method of claim 5 wherein the vector image is a result of Multi -Energy X-ray (micro-, nano-) computed tomography.

7. The method of claim 6 wherein the vector image is a result of Dual -Energy X-ray (micro-, nano-) computed tomography

8. The method of claim 5 wherein the 3D microstructural image is obtained by X-ray fluorescence microtomography.

9. The method of claim 1 wherein the 3D microstructural image is obtained by 3D FIB-SEM.

10. The method of claim 1 wherein the total porosity is measured by gas porosimetry.

11. The method of claim 1 wherein the total porosity is measured by a saturation and weighting technique.

12. The method of claim 1 wherein the total porosity is determined by mercury porosimetry.

13. The method of claim 1 wherein if the sample is polymineral, additionally a 3D mineral mapping of the components is performed.

14. The method of claim 13 wherein the 3D mineral mapping is performed by synchrotron monochromatic X-ray micro- (nano-) tomography.

15. The method of claim 13 wherein the 3D mineral mapping is performed by 3D FIB-SEM with energy dispersive spectroscopy.

16. The method of claim 13 wherein the 3D mineral mapping is performed by combination of 3D X-ray microtomography and 2D scanning electron microscopy with energy dispersive spectroscopy.

17. The method of claim 1 wherein the values in the obtained 3D microstructural image that correspond to the fully void and the fully solid voxels are determined automatically by the computing device.

18. The method of claim 1 wherein the values in the obtained 3D microstructural image that correspond to the fully void and the fully solid voxels are determined manually by an operator.

19. The method of claim 18 wherein the values that correspond the fully void and the fully solid voxels are determined by averaging corresponding group of voxels.

20. The method of claim 18 wherein the values that correspond to the fully void and the fully solid voxels are determined by finding statistical mode in corresponding group of voxels.

21. The method of claim 1 wherein only one of the values in the obtained 3D microstructural image that correspond to the fully void or the fully solid voxels for a particular mineral are determined manually by an operator and other values determined automatically, by the computing device, from condition of equality between experimental total porosity and derived from the overall porosity of the 3D spatial distribution of porosity in the sample.

22. The method of claim Error! Reference source not found, wherein the automatic determination of the values corresponding to the fully void or the fully solid voxels consists in iterative search until satisfaction of the porosity equality condition.

23. The method of claim 1 wherein if the sample is polymineral, physical values for particular minerals required for the normalization and imaged by the 3D microstructural imaging procedure are estimated theoretically.

24. The method of claim 1 wherein if the sample is polymineral, physical values for particular minerals required for the normalization and imaged by the 3D microstructural imaging procedure are measured directly using mineral standards.

25. The method of claim 1 wherein fully solid or fully void standards are artificially added for further determining the values in the obtained 3D microstructural image that correspond the fully void and the fully solid voxels.

26. The method of claim 1 wherein the determined values in the obtained 3D microstructural image that correspond to the fully void and the fully solid voxels are used for creating 3D spatial distribution of porosity in the sample of similar samples imaged under similar conditions.

27. The method of claim 1 wherein the creation of 3D spatial distribution of porosity in the sample includes additional post-processing.

28. The method of claim 277 wherein during the post-processing consists partially solid voxels surrounded by the fully void voxels are replaced with the fully void voxels.

29. The method of claim 277 wherein during the post-processing partially void voxels surrounded by the fully solid voxels are replaced with the fully void voxels.

30. The method of claim 1 wherein the created 3D spatial distribution of porosity in the sample is used for performing numerical simulations of various physical phenomena in the heterogeneous media.