A method and system for improving PORE-scale network modeling of two-phase flow in porous media

Abstract

A method and system for predicting the physical and chemical media properties of a porous media by one or more processing units. The method includes performing a pore-network modeling (PNM) and applying a Mayer-Stowe-Princen (MS-P) semi-analytical theory of drainage to the PNM. The MS-P semi-analytical theory of drainage includes constructing pore boundaries from an image of a pore cross-section. Geometrically and thermodynamically-favorable arc meniscus (AM) are identified from a plurality of AMs. The geometrically and thermodynamically-favorable AM having a contact angle and a curvature radius. An optimal fluid configuration is determined from the geometrically and thermodynamically-favorable AMs.

Description

A METHOD AND SYSTEM FOR IMPROVING PORE-SCALE NETWORKMODELING OF TWO-PHASE FLOW IN POROUS MEDIABACKGROUNDField

[0001] Aspects of the present disclosure generally relate to methods and systems for physical characterization of porous media, and more particularly, to predicting two-phase fluid flow in a disordered porous medium.Description of the Related Art

[0002] Modeling techniques for fluid flow through porous media are broadly implemented for petroleum resource development, materials engineering, food packaging, and medical technology development. Fluid flow modeling techniques may be equipped to illustrate both physical and chemical media properties like permeability, capillary pressure, fluid saturation, contact angle, wettability, or other similar properties, which may be used to characterize fluid behavior within a porous media sample without requiring expensive destruction of the sample.

[0003] Although current techniques for modeling fluid flow through porous media are based on technological advancements made over many years, resultant models may still be tenuous representations of actual porous media. For example, fluid flow models of porous media may require low-resolution implementations to match currently available computational capabilities. As a result, fluid flow models based on porous media having microscale porosities may not accurately reflect physical and chemical properties of the media. Accordingly, there is an impetus to improve the accuracy of fluid flow modeling, including, for example: improving image processing techniques to allow for higher resolution model input and model output, improving image processing techniques to allow for more accurate model input and model output, improving in-situ characterization extraction techniques to better capture fluid behavior in microscale pore features, enhancing computational processing capability to reduce computational expense, enhancing computationalprocessing capability increase modeling speed, increasing automation for iterative modeling steps, improving model capability for dynamic and / or quasi-static modeling of different fluid flow environments, improving model capability for dynamic and / or quasi-static modeling of larger fluid flow environments, and the like.

[0004] Consequently, there exists a need for further improvements in fluid flow modeling of porous media to overcome the aforementioned technical challenges and other challenges not mentioned.SUMMARY

[0005] One aspect of the present disclosure provides a method for predicting two-phase fluid flow in a disordered porous medium by one or more processing units.

[0006] Other aspects provide: an apparatus operable, configured, or otherwise adapted to perform the aforementioned methods as well as those described elsewhere herein; aa non-transitory, computer-readable media comprising instructions that, when executed by a processor of an apparatus, cause the apparatus to perform the aforementioned methods as well as those described elsewhere herein; a computer program product embodied on a computer-readable storage medium comprising code for performing the aforementioned methods as well as those described elsewhere herein; and an apparatus comprising means for performing the aforementioned methods as well as those described elsewhere herein. By way of example, an apparatus may comprise a processing system, a device with a processing system, or processing systems cooperating over one or more message passing interfaces.BRIEF DESCRIPTION OF THE DRAWINGS

[0007] So that the manner in which the above recited features of the present disclosure can be understood in detail, a more particular description of the disclosure, briefly summarized above, may be had by reference to aspects, some of which are illustrated in the appended drawings. It is to be noted, however, thatthe appended drawings illustrate only example aspects and are therefore not to be considered limiting of its scope, may admit to other equally effective aspects.

[0008] Figure 1A is a perspective of an example section of a pore network extracted from porous sandstone, according to embodiments.

[0009] Figure 1B is side view of an example of a high-resolution porous media image, according to embodiments.

[0010] Figure 2 is schematic a flooding instrument, according to embodiments.

[0011] Figure 3 is an example pore network extraction (PNE) process, according to embodiments.

[0012] Figure 4 is a flow chart of a workflow of a MS-P theory based semi- analytical method, according to embodiments.

[0013] Figure 5 is an n-pointed regular star-shaped pore, according to embodiments.

[0014] Figure 6 is a model of the semi-analytical method, according to embodiments.

[0015] Figure 7 is an example of the construction of the arcs that intersect with a pore boundary with a contact angle 0 at both the starting and ending points, according to embodiments.

[0016] Figure 8 is an example of the construction of arcs (dotted curves) that intersect with a pore boundary (solid line) with a contact angle of G at both the starting and ending points, according to embodiments.

[0017] Figure 9 illustrates example available arcs before and after removal of invalid arcs, according to embodiments.

[0018] Figure 10 is a schematic of an Artificial Neural Network (ANN), according to embodiments.

[0019] Figure 11 is a graph of the iteration errors during the estimation of the capillary entry pressure of a pore element using the semi-analytical model, according to embodiments.

[0020] Figure 12 is histogram for the pore-size distributions within of the four rock samples, according to embodiments.

[0021] Figure 13 is a graph of the relative errors between the capillary entry pressures predicted from the semi-analytical method 500 and the corresponding analytical solutions, according to embodiments.

[0022] Figure 14 is a graph of the predicted and measured capillary pressure curves for primary drainage in the Bentheimer sandstone porous media sample, according to embodiments.

[0023] Figure 15 is a graph of the predicted and measured capillary pressure curves for drainage in the Berea sandstone porous media sample, according to embodiments.

[0024] Figure 16 (a)-(c) is the pore fluid-occupancy maps predicted by the enhanced PNM approach for the Bentheimer sample, according to embodiments.

[0025] Figure 16 (d)-(f) is the pore fluid-occupancy maps predicted by the enhanced PNM approach for the Berea sample, according to embodiments.

[0026] Figure 17 is a graph of the predicted capillary pressure curves for MICP test in the Fond du Lac dolostone (FLD) rock sample using the traditional and enhanced pore-network modeling approaches, according to embodiments.

[0027] Figure 18 is a graph of example predicted capillary pressure curves for a MICP test in a Desert Sunset limestone rock sample using the alternative and enhanced pore-network modeling approaches, according to embodiments.

[0028] Figure 19(a) is a visualization of the 3D invading-phase clusters of FLD at a saturation of 0.19, according to embodiments.

[0029] Figure 19(b) is a visualization of the 3D invading-phase clusters of FLD at a saturation of 0.41, according to embodiments.

[0030] Figure 19(c) is a visualization of the 3D invading-phase clusters of FLD at a saturation of 0.61, according to embodiments.

[0031] Figure 19(d) is a visualization of the 2D fluid-occupancy maps at a saturation of 0.19, according to embodiments.

[0032] Figure 19(e) is a visualization of the 2D fluid-occupancy maps of FLD at a saturation of 0.41 , according to embodiments.

[0033] Figure 19(f) is a visualization of the 2D fluid-occupancy maps of FLD at a saturation of 0.61 , according to embodiments.

[0034] Figure 20(a) is a visualization of the 3D invading-phase clusters of DSL at a saturation of 0.21, according to embodiments.

[0035] Figure 20(b) is a visualization of the 3D invading-phase clusters of DSL at a saturation of 0.42, according to embodiments.

[0036] Figure 20(c) is a visualization of the 3D invading-phase clusters of DSL at a saturation of 0.61, according to embodiments.

[0037] Figure 20(d) is a visualization of the 2D fluid-occupancy maps of DSL at a saturation of 0.21 , according to embodiments.

[0038] Figure 20(e) is a visualization of the 2D fluid-occupancy maps of DSL at a saturation of 0.42, according to embodiments.

[0039] Figure 20(f) is a visualization of the 2D fluid-occupancy maps of DSL at a saturation of 0.61 , according to embodiments.

[0040] Figure 21 is a graph of the Pearson correlation between the geometric features and the capillary entry pressure of pores, according to embodiments.

[0041] Figure 22 is a scatterplot matrix of the geometric features against the capillary entry pressure, according to embodiments.

[0042] Figure 23 is a graph between the predicted capillary entry pressures using the ANN and the corresponding capillary entry pressures, according to embodiments.

[0043] Figure 24 is an example device for characterizing pore networks, according to embodiments.

[0044] To facilitate understanding, identical reference numerals have been used, where possible, to designate identical elements that are common to the figures. It is contemplated that elements and features of one aspect may be beneficially incorporated in other aspects without further recitation.DETAILED DESCRIPTION

[0045] In the following, reference is made to aspects of the disclosure. However, it should be understood that the disclosure is not limited to specifically described aspects. Instead, any combination of the following features and elements, whether related to different aspects or not, is contemplated to implementand practice the disclosure. Furthermore, although aspects of the disclosure may achieve advantages over other possible solutions and / or over the prior art, whether or not a particular advantage is achieved by a given aspect is not limiting of the disclosure. Thus, the following aspects, features, embodiments, and advantages are merely illustrative and are not considered elements or limitations of the appended claims except where explicitly recited in a claim(s). Likewise, a reference to “the disclosure” shall not be construed as a generalization of any inventive subject matter disclosed herein and shall not be considered to be an element or limitation of the appended claims except where explicitly recited in a claim(s).

[0046] TThhee present disclosure rreellaatteess to techniques ffoorr physical characterization of porous media. Specifically, the techniques discussed herein may be implemented for non-destructive extraction of geometry information, capillary pressure and wettability information from porous media samples. The porous media sample may comprise a digital rock sample, a rock sample, a core sample, a plastic sample, a tissue sample, a disordered space, a rough-walled fracture space, or any other organic or inorganic sample having pore space ascertainable through imaging techniques.

[0047] A thorough grasp of fluid flowthrough porous spaces of certain materials may be consequential to enhancing technical efficacy of fluid flow techniques in a wide range of industries. Models of fluid flow are useful to describe physical and chemical characteristics of a porous material and may help to highlight the material’s optimal usage. Often, networks of pores within a material are extremely small (e.g., microscale in size). Techniques for characterizing these pore networks are hindered by the computational expense of modeling at a microscale. To alleviate computational burdens, pore network modelling techniques often implement generalized characterization techniques at expense of model accuracy. Extrapolation errors caused by such imprecise characterization may result in mischaracterization of physical and chemical characteristics of the porous material. In many cases, these errors render such models impractical for regular use. Accordingly, ideal modeling of fluid flow through porous media would allowfor rapid, accurate characterization of microscale pore spaces that may be performed without inhibitive computational expense.Introduction to Pore Network Modeling

[0048] Modeling techniques for fluid flow through porous media may illustrate both physical and chemical porous media properties. Models of porous media may be used to ascertain permeability, capillary pressure, fluid saturation, wettability, buoyancy, and the like to a greater degree of accuracy more comparable to physical flooding of a porous media sample. Additionally, physical and chemical properties determined using pore network modeling techniques may be used to characterize in-situ fluid behavior as it travels through the porous media under a wide variety of wettability and flooding conditions. These conditions are not typically accessible to uusseerrss performing conventional physical flooding characterization techniques.

[0049] Permeability is the tendency of the porous media to allow liquids to flow through the porous media. Capillary pressure is the pressure difference existing across the interface separating two immiscible fluids. Fluid saturation is the measurement of fluid present in the pore spaces of the porous media. Contact angle is a measured angle between a fluid-fluid or a fluid-gas interface at a point where it meets a solid surface. Wettability is the ability of a liquid to maintain contact with a solid surface. Wettability may vary depending on wettability conditions and the type of wetting liquid present in the porous media sample. For example, a water-wet medium may show a lower wetting affinity to the oil phase than an oil-wet medium, where higher or lower wetting is determined with respect to a given phase. In certain cases, the correlation between wettability and viscosity ratio may not be straightforward, as there may be water or oil wet conditions with similar viscosities.

[0050] Figure 1A is a perspective of an example section of a pore network extracted from porous sandstone. A modeled pore network is a practical description of a porous medium targeted for fluid flow modeling. The section of the pore network describes the porosities of various size and shape present in that portion of the sandstone, and may be used to model fluid flow through thoseporosities for various wettability conditions. Three-dimensional (3D) portions of a pore network model may more accurately characterize the porous media sample either alone or in combination with other 3D portions of the pore network model.

[0051] Pore network models (e.g., of Figure 1A) may be extracted from images of a targeted porous medium and used to model multi-phase fluid flow using physically-based displacement mechanisms (PBDMs) across pores defined in a pore network. PBDMs may represent an estimated displacement of a modeled fluid in response to movement of another fluid or gas within the pore network. As immiscible phases react with one another throughout the pore network during fluid flooding, PBDMs are induced where, for example, capillary pressure across a meniscus exceeds the wettability constraints on either phase. Fluid saturation, contact angle, buoyancy, and the like may also affect PBDMs throughout a pore network. By utilizing a pore network model extracted from a porous media sample, a user is able to ascertain PBDMs through the porous media sample under a wide variety of wettability conditions in order to ultimately obtain, for example, useful permeabilities for a larger sample of the porous medium without degrading a porous media sample via repeated physical flooding.

[0052] Figure 1B is side view of an example of a high-resolution porous media image. The high-resolution porous media image is taken by a scanning instrument from a single rock sample and segmented for characterization. To properly generate PBDMs at a pore scale for the targeted porous media, imaging may capture complex geometries of the targeted porous media at a resolution sufficiently high to retain acceptable accuracy. Pores may be defined as a complex polyhedron having at least a center 102 and spherical and effective diameters. Connective throats 104 between pores may also be defined. In many cases, image resolution may be in micrometers to capture applicable pore detail. High-resolution pore models enable an accurate rendering of the fluid flow characteristics described above as ascertained at each pore and for each PBDM.

[0053] PBDMs may occur upon flooding or draining of a pore network model, where aqueous phase injection or removal is iteratively simulated through the pore network. Aqueous flooding and aqueous draining may be implemented in variousmodeled wettability conditions, where certain fluids are present prior to the start of a simulation. Wettability conditions may include at least water-wet, oil-wet, or mixed-wet conditions. During aqueous flooding, injected water may displace immiscible fluid preexisting in the pore network model. During aqueous draining, injected immiscible fluid may displace water preexisting in the pore network model. In certain cases, flooding and draining may be fluid flooding and fluid draining. In some cases, the fluid is oil.

[0054] Figure 2 is schematic a flooding instrument 200. Flooding or draining of a pore network model may be simulated based in part on scanned images of physical flooding implemented by the flooding instrument 200. In some cases, a porous media may undergo a core-flooding experiment to establish an irreducible water saturation, a residual oil saturation, or both. Core-flooding is enabled by a set of pumps 202, rupture disks 204, pump lines 206, 208, 210, 212, and 214, differential pressure transducers 216, and source buckets 218. 220, and 222 working in tandem to flood a porous media sample loaded in a core holder. In some cases, a scanning instrument (e.g., a micro computed tomography (micro-CT) scanner) captures a dry reference image prior to flooding. Scanning occurs in a field of view defined within the core holder. The porous media sample may be flooded with brine from bucket 220 via the brine tubing line 206 and scanned again to ensure that the porous media sample is fully saturated. Once the brine flooding is complete, the absolute permeability of the porous media sample may be obtained. The oil flooding may be performed alongside additional brine flooding. Any fluid expelled as a result of overburden pressure (i.e., pressure that compacts pore space and reduces permeability) is transported via the confining fluid line 208 and collected in bucket 222. Any fluid expelled as a result of the flooding procedure is transported via the effluent fluid line 212 and collected in bucket 224. In many cases, core sample pressure may be iteratively adjusted during flooding. Pressure may be recorded by one or more differential pressure transducers 216 coupled to the core holder via a transducer line 214.

[0055] Scanned innages obtained from flooding procedures performed by the flooding instrument 200 may be used to extract a pore map representative of theporous media sample prior to, during, or after flooding. The images may be processed to determine characteristics of fluid flow through the porous media sample. In many cases, the images may also be used to extract a representative pore network model.

[0056] Imaging of porous media is typically performed using micro-CT imaging. In many cases, commercial micro-CT scanners (e.g., Zeiss scanners) are available for imaging necessary to perform pore network modelling. Images of porous media taken by micro-CT scanners are at a sufficiently high resolution to create a microscale digital image of the porous media.

[0057] In the current state of the art, there exists a challenge of extracting porous media characteristics in a manner precise and repeatable to ensure the ultimate stability of future simulations. Currently, techniques for porous media characterization require lengthy step-wise processing known to incur undue computational expense and introduce instability to characterize of the porous media sample. As a result, users may not be able to rely on characterization output to simulate flow conditions in a useful way.Aspects Related to Pore-Network Modeling of Two-Phase Flow in Mixed-Wet Media

[0058] Fluid flow modelling through porous media is often utilized to enhance petroleum resource development. In recent years, global demand for energy resources has mobilized development of unconventional petroleum reservoirs as targets for hydrocarbon extraction. Example geological formations that comprise these newly developed hydrocarbon reservoirs are ultra-tight shale formations resistant to primary petroleum extraction techniques. A matrix of an ultra-tight unconventional shale reservoir may be characterized by low permeability and low porosity. To extract hydrocarbons from the ultra-tight shale matrix, secondary and tertiary petroleum extraction techniques seek to maximize oil production through the microscale pore networks that comprise a substantial amount of the porosity in the shale matrix.

[0059] A robust understanding of fluid flow through microscale pore networks of hydrocarbon reservoirs may be consequential to extracting the trillions of barrels of oil and gas still housed in shale formations globally. Models of fluid flow througha pore network that incorporates permeability, capillary pressure, fluid saturation, contact angle, and wettability may help to elucidate specific steps to be taken during resource development to optimize petroleum production. Even so, techniques for characterizing these microscale pore networks are hindered by the computational expense of modeling microscale pore networks and extrapolation errors caused by oversimplified characterization of pore geometries.

[0060] As discussed above, ideal modeling of fluid flow through porous media would allow for precise, quick, and repeatable characterization of a porous media sample. In a case where the porous media sample is, for example, a cylindrical core sample of a rock having a length of six inches and a diameter of one inch, the core sample is likely to have porosity and permeability characteristics that vary across its length and width. This is common in core samples and especially in core samples representative of ultra-tight oil formations. Geological processes that form certain oil-bearing rocks can produce heterogeneous morphological features in the rock that may be present even at a micrometer scale. This is especially true for oil- bearing carbonate rocks, which contain micro-porosities that contribute significantly to the overall porosity of the rock. These microscale morphological features may affect the pore network of the core sample, altering the porosity and permeability throughout a core sample. Thus, accurate characterization of fluid flow through a core sample may depend on precisely ascertained and verifiable microscale geometries sufficient to detect heterogeneous properties of a pore network. Using conventional estimation techniques that cannot consistently capture the heterogeneity and complexity of either the core sample or the fluid- fluid interfaces present therein may result in characterization of a porous media sample that cannot be used to consistently describe fluid flow through the core sample.

[0061] According to aspects of the present disclosure, the pore network construction procedure described herein may be performed by a processing system architecture comprising at least one or more central processing units (CPUs) operating independently or in combination with one or more graphics processing units (GPUs). Techniques for pore network construction performed byprocessors operating in sequence may be unable to support the computational load associated with microscale pore networks. Techniques implemented herein may address issues with resolution sufficiency and issues arising from the computational expense of high-resolution pore network extraction. The one or more CPUs and / or the one or more GPUs may perform the procedures according to a non-transitory computer readable medium that causes the one or more CPUs and / or the one or more GPUs to perform any and all steps of the procedure. Each of the one or more CPUs may be utilized in combination with a memory having the computer readable medium stored thereon. Each of the one or more CPUs may be utilized in combination with one or more processors. Each of the one or more processors may be parallel processors. Each of the one or more GPUs may be utilized in combination with a memory having the computer readable medium stored thereon. Each of the one or more GPUs may be utilized in combination with one or more processors. Each of the one or more processors may be parallel processors. Each of the CPUs and the GPUs may operate independently, or may operate using a message passing interface (MPI) enabling communication between one or more parallel processors for performing the procedure. This may include CPU-CPU communication, CPU-GPU communication, and / or GPU-GPU communication.

[0062] Pore-network modeling (PNM) has served as a practical simulation tool for predicting multi-phase flow functions of porous media, owing to its superior computational efficiency. However, in traditional pore-network models, the rock pore space is simplified, and the fluid flow and solute transport are solved in a network of pore elements with idealized geometries. Such simplification can lead to inaccurate predictions when the original pore space features highly diverse and complex geometries as those observed in carbonates. To address this issue, aspects described herein propose to couple the PNM approach with a new computationally efficient semi-analytical method that can faithfully predict capillary entry pressure and the corresponding fluid configuration of piston-like displacements using real, two-dimensional (2D) images of the pore space. To this end, a procedure for identifying and extracting 2D pore cross-sections from thecorresponding 3D images of the pore space is developed and incorporated within an existing pore-network extraction platform. Aspects of the present disclosure provide means for analysis of geometric features of the pore space and their effects on the capillary entry pressure of pores. For example, critical geometric features of hydraulic radius, circumradius, sphericity, and area have strong associations with the pore capillary entry pressure. Aspects of the present disclosure provide an Artificial Neural Network (ANN) to predict the capillary entry pressure of pores using their critical geometric features. In some cases, one or more processors may sample an ensemble of about 80,000 pores to about 120,000 pores (e.g., 100,000 pores) to characterize their geometric features, and compute their capillary entry pressures using the semi-analytical method 500 to train the ANN. The ANN procedures show excellent predictive accuracy (e.g., accuracy having an R2of about 0.995).

[0063] Aspects of the present disclosure provide an efficient and robust framework to carry out pore-scale displacements based on the fluid pressure fields that are updated frequently during flow simulations. The efficiency of framework combined with the parallelization of the platform across one or more CPUs may facilitate performing fluid flow prediction in large core-sized pore networks within a practical amount of time.

[0064] With the rapid development of imaging techniques and scientific computations, Digital Rock Technology (DRT) has become a powerful tool for investigating pore-scale physics and predicting the properties of multi-phase flow in porous media involved in many subsurface applications. In general, pore-scale flow models of DRT can be categorized into Direct Numerical Simulations (DNS), and Pore-network Modeling (PNM). In DNS, the governing equations of fluid flow and solute transport processes are solved directly on a discretized pore space. Due to the complexity of pore space boundaries and pore-scale fluid dynamics, the computational cost of this approach is prohibitively expensive. As an alternative, the PNM approach is developed by representing the pore space as a network of pores and throats with idealized geometries, which significantlyimproves its computational efficiency, thereby extending its application to much larger porous media than those studied by DNS.

[0065] In general, PNM starts with constructing a representative pore network from high-resolution images of the porous media. By design, pore networks can capture the comprehensive topological configurations and connectivity details of the pore space. The PNM approach allows prediction of the flow and transport properties (e.g., capillary pressure and dispersion coefficient) of porous media samples. However, its application in more complex porous materials (e.g., carbonates) is still limited. This is mainly because the simplified representation of the complex pore space as basic convex polygons (circles, triangles, and rectangles) can reduce the accuracy of the estimated pore-level properties that govern multiphase flow in such porous media.

[0066] The development of coupled pore-scale models is not yet practical for a direct pore-by-pore analysis in large-sized porous media due to the computational cost associated with the underlying DNS. To this end, a thermodynamically consistent method known as the Mayer-Stowe-Princen (MS-P) semi-analytical theory of drainage may be leveraged to accurately and efficiently predict pore-level two-phase flow properties. Based on a rigorous energy-balance analysis, the semi- analytical MS-P method has a distinct computational advantage over DNS.

[0067] To better account for the effect of pore geometries on capillary entry pressures at the pore level, a workflow couples the PNM with the MS-P method. Implementations first extend an existing seamless PNE platform via one or more processors to identify and extract the 2D cross-sections of the pore spaces from their high-resolution micro-CT images. Next, implementations described herein introduce an improved MS-P platform for modeling two-phase displacement under varying wettability conditions and its parallel implementation. Subsequently, implementations described herein couple the MS-P platform with PNM to predict capillary pressure curves and fluid-occupancy maps of two-phase displacement. Examples provided illustrate such implementations in both naturally occurring sandstones and carbonates using real pore geometries. Finally, implementationsdescribed herein utilize ANN such that the pore capillary entry pressures may be quickly predicted with known critical geometric features.

[0068] Figure 3 is an example pore network extraction (PNE) process. The example PNE process is performed by one or more processors according. A three- dimensional micro-CT image is decomposed into overlapping images (left). Distinct networks are extracted in parallel, using medial axis analysis and graph reduction techniques to guarantee consistency across the region of overlap (middle). A complete pore network is constructed using a seamless merging technique to combine the extracted subnetworks (right). In the ball and stick representation of each pore network, purple spheres and yellow cylinders are used to represent pore bodies and throats, respectively.

[0069] At the start of the seamless PNE process, the global image is decomposed into a set of overlapping sub-images via one or more processors, from which a set of corresponding sub-networks are extracted in parallel. Afterward, these sub-networks are seamlessly merged by one or more processors to obtain a cohesive representative description of the imaged porous media. The processors then implement a graph and medial-axis-driven PNE method to generate each sub-network. This approach may optimize the consistency of these sub-networks over each region of overlap, allowing adjacent sub-networks to be seamlessly merged. A distinct advantage of this approach is that the voxel information of pore elements is preserved, enabling the direct determination of their cross-sections and other morphological properties.

[0070] In some cases, the cross-section of a throat element is defined along the mid-plane that passes through its center and is perpendicular to the line connecting the centers of its two neighboring pore bodies. As an example, for a throat element, with center xc, connected to pore bodies i and j with respective centers xiand xj, the mid-plane is defined as the set of points x that satisfy the constraint equation shown in Equation 1 :

[0071] The cross-section of a pore body is constructed similarly, except the normal vector of the mid-plane is defined as the resultant of vectors pointing fromthe center of the pore body to its adjacent throats. That is, for a pore body i, with center xi, we define its corresponding mid-plane as the set of points x satisfying the constra

[0072] where Ω1represents the set of pore throats attached to i.

[0073] After the mid-plane has been defined within a search radius from the center of the pore element, the one or more processors identify all the voxels of the segmented micro-CT image that intersect with the mid-plane. To ensure that the cross-section captures all the voxels within the boundaries of the pore element, the processors use a search radius that is a multiple of its spherical radius. The set of identified voxels in 3D space may be mapped to the mid-plane to construct the 2D cross-sectional image of the pore element. The resolution of the resulting image may be improved by refining the discretization that is traversed during this mapping process. Finally, artifacts within the cross-sectional image may be removed using a sequence of morphological image processing steps, including image dilatation and erosion operations and hole-filling procedures. The computed cross-section may be then output as a binary 2D image for use in the enhanced PNM.

[0074] In some cases, an improved implementation of MS-P analysis is implemented to compute the capillary entry pressure of piston-like displacements taking place in pore elements. Equation 3-9, as discussed below, show the critical foundation of the MS-P approach. Based on an energy-balance description, the generalized MS-P analysis states that for an infinitesimal two-phase movement in a closed system at a constant temperature and a fixed total volume, the change in Helmholtz free energy of the system, dF, is the summation of the work done on the system and the change in its interfacial free energy, shown in Equation 3:

[0075] wrespectively; pi and Vi are the pressure and volume of the bulk phase i, and crij and Aij are the interfacial tension and interfacial area between phases i and j, respectively. For a reversible displacement process at constant temperature andvolume when the system is at equilibrium with a thermodynamically favorable fluid configuration, its Helmholtz free energy is minimum, i.e., dF = 0, is shown in Equation 4:

[0076] where Pc is the pressure difference between the wetting and non-wetting phases.

[0077] The interfacial tension between the two fluids may be defined by Equation 5:

[0078] where 6 is the contact angle between the fluid interface and the solid surface. By combining Equation 4 and Equation 5, the capillary pressure equation may be formulated as Equation 6:

[0079] If the pore space is assumed to have a constant cross-sectional shape, for a movement dx of the invading fluid, it may be inferred that Equation 7 applies:

[0080] where Lnwand Lnsare the total contact lengths of the non-wetting- wetting and non-wetting-solid interfaces, respectively, at the thermodynamic equilibrium condition. Then, Equation 6 can be further simplified as Equation 8:

[0081] Provided that one or more arc menisci (AMs) form at the corners of the pore element, their radius, r, can be characterized using at least Equation 8, as Equation 9:

[0082] As suggested by Equation 9, the capillary entry pressure of the piston- like displacement in a pore element with any geometric shape can be analytically computed if the wettability condition (0), the fluids interfacial tension (a), and the geometric variables corresponding to the fluid configuration areknown or ascertainable. However, except for certain regular-shaped (e.g., triangle, square, and star) pores, the optimal fluid configuration at the capillary entry pressure cannot be determined a priori. In such cases, at the condition for a piston- like invasion, the drainage curvature radius of the invading interface in a pore element may be related to the effective geometric parameters of the pore space and the fluid configuration according to Equation 9.

[0083] Figure 5 is an n-pointed regular star-shaped pore. The total contact lengths of the non-wetting-wetting and non-wetting-solid interfaces, Lnwand Lns, are given by Equation 10 and Equation 11:

[0084] where 0 is the contact angle; β is the half corner angle; and P is the perimeter of the star.

[0085] The effective area that is occupied by the invading phase can be obtained through a geometric analysis:

[0086] And the total area of the star, A, equals:

[0087] where rins is the inscribed radius of the star.

[0088] By combining Equations 9-14, and letting

[0089] After re-arranging the terms, a quadratic function of the curvature radius of the invading interface is expressed as Equation 15:

[0090] Note that

[0091] Finally, the curvature radius, r, is characterized as the valid root ofEquation 15, shown in Equation 16:

[0092] Figure 4 is a flow chart of a workflow of a MS-P theory based semi- analytical method 500. Figure 6 is a model of the semi-analytical method 500. For pore elements with complex, irregular cross-sections, Equations 3-9 of the MS-P cannot be solved directly in an analytical manner, therefore, the semi-analytical method 500 is utilized. The semi-analytical method 500 for predicting capillary entry pressure of pores contains three major operations which may be performed by one or more processors: operation 410 - construction of the pore boundary from the image of the pore cross-section; operation 440 - identification of geometrically-feasible and thermodynamically-favorable AMs from possible arcs obtained with the given contact angle and curvature radius; and operation 450 - determination of the optimal fluid configuration from the identified AMs.

[0093] Operation 410 (pore boundary construction) includes operation 412, 414, and 416. At operation 412, as shown in (a) of Figure 6, a high-resolution image of the pore cross-section image is obtained. At operation 414, the pore boundary is extracted from the high-resolution image. One or more processors use a computer vision and image analysis library (e.g., OpenCV-Python) to extract the boundary points of the pore space from its segmented cross-sectional image. At operation 416, as shown in (b) of Figure 6, the extracted pore-boundary is smoothed using a moving-average filter. The processors apply the moving- average filter along the extracted boundary points to smooth out artificial discontinuities of the boundary curve arising from the pixelated image. The set of the extracted boundary points as B.

[0094] At operation 420, as shown in (c) of Figure 6, an initial guess of the curvature of the arc meniscus (AM) is performed. Figure 7 is an example of the construction of the arcs that intersect with a pore boundary with a contact angle 0 at both the starting and ending points.

[0095] The initial guess includes all the arcs intersecting with the pore boundary with contact angle 9. The iterative procedure starts with an initial guess of the AM’s radius that is computed by assuming the invading fluid will occupy the entirepore space after the piston-like invasion. In this regard, Anw and Lns in Equation 9 are substituted by Apotion 17:

[0096] wwhheerree Aporeand Lporemay be the cross-sectional area and perimeter length of the pore cross-section.

[0097] The one or more processors may then generate a set of all the AMs with radius equal to rinitthat form the contact angle 9 at their two intersecting points with the pore boundaries. To do so, first, a circle with a radius of rinit may be traversed along the pore boundary points (B) in two opposite directions, for example, clockwise and counterclockwise directions, to locate all possible arcs that can form the contact angle 9 at either their starting or ending intersecting points. The center of these arcs can be determined through the following geometric relations:

[0098] where xi and yi may be the coordinates of pore boundary at point i. The orientation angle, a, of the boundary point i may be the angle measured between the x-axis and the boundary tangent line at i in the counterclockwise direction. By finding the intersection points of two closed curves, the one or more processors may obtain the centers of all the arcs forming the contact angle 9 at both ends. The centers of the arcs forming a valid contact angle at their starting and ending points are shown as lighter dots and darker dots, respectively, while the centers of arcs forming valid contact angles at both ends are shown as thick, dark dots.

[0099] At operation 430, as shown in (d) of Figure 6, geometrically- and physically-allowed arcs are determined. Operation 440 (valid arc meniscus (AM) identification) includes operations 441-445. In some cases, many of the arcs may not form physically meaningful AMs. Thus, a key component of the method 500 identifies and eliminates the nonphysical arcs that cannot form stable AMs. In general, there are four types of invalid arcs, including ones that i) are located outside the pore boundary, ii) have diverging curvatures and thus are physicallyunstable, iii) intersect with existing arcs, and iv) coexist with other arcs in the same corner or neck region of the pore space.

[0100] At operation 441 , arcs not forming a valid contact angle are removed. At operation 442, arcs having diverging curvatures and, thus, are physically unstable, are removed. At operation 443, arcs outside of the pore boundary are removed. At operation 444, extra arcs intersecting with other arcs are removed. At operation 445, the AM types are identified and the possible fluid configurations are illustrated.

[0101] To address uncertainties inherent to the invalid arc types, the one or more processors may implement thermodynamics-based criterion to identify the most stable arc through an energy balance analysis. Figure 8 is an example of the construction of arcs (dotted curves) that intersect with a pore boundary (solid line) with a contact angle of θ at both the starting and ending points.

[0102] In one example, for a pore corner depicted in Figure 8, arcs A1A2, B1B2, and C1 C2 may remain after preliminary analyses are conducted to remove other invalid arcs of types i-iii. In other words, the drainage displacement may leave behind a wetting corner in a closed region defined by A1OA2, B1OB2, or C1OC2 together with their corresponding arcs. Assuming that AF may be extra energy allowing elimination of the entire wetting corner after the displacement, according to Equations 3 and 7, AF can be characterized as shown in Equation 19:

[0103] Wheremay equal the corner area in front of the arc Lc^nermay equal the total contact length between the wetting phase and the solir may be the radius of the arc. The most stable arc will correspond to the one that may allow the largest AF to be eliminated. Therefore, this corner arc may be preserved while the rest of the arcs existing in the same corner area may be removed.

[0104] At operation 446, extra arcs located in the same pore corner or neck region are removed. Figure 9 illustrates the available arcs before and after operation 446. To test the proposed thermodynamics-based criterion, primary drainage simulations in a selected pore space, as shown in Figure 9, under two different contact angles of (e.g., 20° and 70°) representing the strongly water-wet and weakly water-wet conditions, respectively. As shown in (a) of Figure 9 for the strongly water-wet case (0 = 20°), many arcs are left at the corners of the pore space after invalid arcs of types i-iii are eliminated. The abundance of the remaining arcs is due to the strong adhesive force between the wetting phase and the rock surface, which in return causes the invading non-wetting phase to minimize the region it can occupy. As shown in (b) of Figure 9, by applying this criterion to remove the extra arcs coexisting with others in the same corner region, only the innermost ones that minimize the total Helmholtz free energy of the system are preserved. This result is consistent with the previous computational observations that are obtained under the condition of 9 = 0. As shown in (c) of Figure 6, in the weakly water-wet case, where the contact angle is 70°, the potential for arc formation at the pore corners decreases as the rock surface becomes more neutral wet towards the invading phase. As shown in (b) of Figure 9, when identifying invalid arcs of type iv, the geometric approach should retain the innermost arc A1 A2, even though under such wettability conditions, preserving arc B1 B2 can result in a more stable fluid configuration based on the energy balance analysis.

[0105] Depending on their locations, the remaining arcs may be classified into two types of AMs, corner and neck AMs. As illustrated in (d) of Figure 9, the corner AMs are the arcs behind which wetting corners remain after the piston-like displacement. Here, we define the back of an arc as the side that does not contain the arc center. The neck AMs are pairs of arcs that form in the same pore constriction region but are in opposite directions. A pair of neck AMs, AMi and AMj, may satisfy Equation 20:

[0106] Where BAMi and BAMj may represent pore boundary points in the front of the AMi and AMj, respectively, and Bpore may represent the full set of the pore boundary points.

[0107] Operation 450 (optimal fluid configuration determination) includes operation 452, operation 454, and operation 456. At operation 452, as shown in (e) of Figure 6, the phase regions are identified. The AM configurations include two oil regions O1 and 02 and one water neck region W1 are identified, resulting in four possible piston-like displacement configurations. Depending on the required capillary entry pressure, the invading fluid may occupy region O1 , 02, O1 + 02, or O1 + W1 + 02 that are enclosed by different combinations of the AMs and pore boundary segments. For each of the possible fluid configurations, the corresponding capillary entry pressure is computed according to Equation 8 and the one with the minimum capillary entry pressure is determined as the most optimal fluid configuration.

[0108] At operation 454, as shown in (f) of Figure 6, the optimal fluid configuration of the displacement is identified. The semi-analytical method 500 may be implemented by one or more processors to determine the optimal fluid configuration after displacement. The different possible invading and defending phase regions are first constructed by traversing through all the pore boundary points separated by the two types of AMs. At operation 456, the AM curvature is updated.

[0109] At operation 460, a query regarding whether the AM curvature has converged is performed. If the AM curvature has converged, the semi-analytical method 500 writes out the predicted capillary entry pressure and optimal fluid configuration and the semi-analytical method 500 ends. If the AM curvature has not converged, the method 400 proceeds to the next iteration starting from operation 430 to identify new valid AM(s) with the updated capillary pressure and AM's curvature radius. The method 500 and the corresponding optimal fluid configuration repeat until the curvature radius converges according to Equation 21 :

[0110] Where ε may be a relative error tolerance and i may be the iteration number. In general, the semi-analytical method 500 may produce a model that converges rapidly; for example, where an error tolerance of about 10-7is satisfied after about 5 to about 10 iterations (e.g., about 7 iterations) for the pore shown in Figure 11 , and the entire computational procedure may about 1 second to about 10 seconds (e.g., about 3.9 seconds) on a single thread with a processing speed of about 1GHz to about 5GHz (e.g., about 2.2 GHz). In some cases, the model converges faster for weaker water-wet conditions because a fewer number of geometrically-allowed arcs can be initially formed. Overall, the semi-analytical method 500 produced using aspects of the present disclosure is computationally efficient.

[0111] Figure 11 is a graph of the iteration errors during the estimation of the capillary entry pressure of a pore element using the semi-analytical model. The pore cross-section is extracted from a Bentheimer sandstone sample and the input contact angle is 20°. Examples of the evolving fluid configurations are shown for selected iterations.

[0112] Figure 12 is histogram for the pore-size distributions within of the four rock samples. In some cases, the pore networks extracted from porous media samples according to the method 500 described herein may be especially large. Specifically, such pore networks may have a substantial numbers of pore elements (e.g., as large as about 150,000 pore elements to about 250,000 pore elements, though other values are contemplated). Consequently, the computational workload for predicting capillary entry pressures of all the pore elements may be significant. To reduce the computational time, the semi-analytical method 500 may be executed on processors running in parallel using a high-level interface (e.g., the Python interface, concurrent.futures), designed for asynchronously executing callables.

[0113] According to certain aspects of the present disclosure, the one or more processors may implement a PNM approach to simulate the two-phase primary drainage displacements in the porous media. Specifically, all the extracted pore networks of the porous media may be initially filled with the wetting fluid. For theprimary drainage process, the non-wetting fluid may be introduced at the inlet of the network under a constant injection pressure, while both fluids can be produced at the outlet of the network with a constant production pressure. The injection pressure of the invading fluid may be gradually increased until either a target saturation is achieved or no more pore-scale displacement can take place. The invasion processes may be assumed to be quasi-static, meaning that the viscous pressure drop of the flow is negligible. A no-flow condition may be applied on the side boundaries of the pore networks.

[0114] Primary drainage simulations may be conducted by one or more processors in both enhanced and alternative pore networks of the rock samples. In the enhanced PNM approach, the semi-analytical workflow may first be applied to all the pore and throat elements to predict their capillary entry pressures using the real pore and throat cross-sections. In an alternative PNM approach, the capillary entry pressures of pores and throats are estimated using the analytical solutions derived for the simplified pore geometries.

[0115] Table 1 is the definition of features. In addition to the inscribed radius and shape factor that are typically adopted in pore-network models for estimating capillary entry pressures of pores, the one or more processors may collect twelve more geometric features and investigate their relations with the capillary entry pressures of pores having arbitrary geometries, according to certain aspects of the present disclosure. These features include area, aspect ratio, circumradius, convexity, equivalent diameter, extent, hydraulic radius, mean distance, perimeter, roundness, solidity, and sphericity.Table 1. Definitions of Features.

[0116] Although the computational efficiency of the semi-analytical method 500 herein may be much higher than DNS, running the semi-analytical method 500 on increasingly large pore networks may still be expensive. In addition, the general relationship between the capillary entry pressure and the pore geometric features is complex and nonlinear. In view of these, the one or more processors may construct and train a sequential artificial neural network (ANN) to predict the capillary entry pressure of pores using their geometric features.

[0117] Figure 10 is an Artificial Neural Network (ANN) 1000. The ANN 1000 architecture includes an input layer, hidden layer(s), an output layer with a single neuron, and a Batch Normalization (BN) layer that may be added in front of every hidden layer to enhance the converging stability of the network. ANN 1000 contains numerous adjustable hyperparameters, such as the learning rate, number of hidden layers, and number of neurons in each layer, all of which have large impacts on the model accuracy. To effectively explore the hyperparameters space, the processors may apply randomized search cross-validation to obtain the optimal combination that yields the best performance on a validation data set.

[0118] The ANN 1000 may be constructed, trained, and evaluated using a high- level deep learning methodology. The one or more processors may first split the entire data set into the training, validation, and test subsets, all of which are normalized by applying the Min-Max scaling. With a thorough assessment of the performances of the ANN 1000 with several types of optimization strategies, the one or more processors may implement a Sigmoid function as the activation function, the He initialization as the weight initializer, and the Adam solver with Nesterov momentum as the optimizer. To prevent overfitting, an early-stopping technique may be adopted to terminate the training process once the validation error reaches a minimum. During the training process, the processors may compute the mean squared error (MSE) loss function to monitor the errors between the labels and the corresponding predictions.

[0119] While the MS-P theory based semi-analytical method 500 provides a way of estimating capillary entry pressure of pores directly using cross-section images, the ANN 1000 can compute the capillary entry pressures of a large amount of pores with various geometric properties. This enables the training of a machine learning (ML) based ANN 1000. With a well-trained ANN 1000, the capillary pressure of any pore with known geometric properties can be estimated. In other words, the MS-P theory based semi-analytical method 500 can enable the generation of data to train the ANN 1000, which can then estimate capillary entry pressure of any pore without explicitly going through the method 500.

[0120] Figure 13 is a graph of the relative errors between the capillary entry pressures predicted from the semi-analytical method 500 and the corresponding analytical solutions. Three idealized pore conditions under various wettability conditions are used to determine the relative errors: equilateral triangle; square; and five-pointed star with a corner angle of 36°. The semi-analytical method 500 may output a fluid configuration corresponding to fluid flow through a target porous media sample. To obtain and analyze such output, the one or more processors may first conduct numerical simulations of oil-water piston-like displacements in idealized pore elements with the same inscribed diameter but different cross- sectional geometries.

[0121] During primary drainage displacement in these angular pore geometries, wetting layers may form in the corners of the pore space. Accordingly, the semi-analytical method 500 reproduces the expected fluid configurations with wetting phase residing in the corners of these pores. The processors may also compute the relative errors between the predicted capillary entry pressures using the semi-analytical method 500 and the corresponding analytical solutions under different wettability conditions. The relative errors are in the range of about 10"5to about 102, indicating the accuracy of the semi-analytical method 500 over a wide range of contact angles.

[0122] Figure 14 is a graph of the predicted and measured capillary pressure curves for primary drainage in the Bentheimer sandstone porous media sample. Figure 16 (a)-(c) is the pore fluid-occupancy maps predicted by the enhanced PNM approach for the Bentheimer sample. The one or more processors perform primary drainage simulations in Bentheimer sandstones and compare the predicted capillary pressure curves with the experimental data corresponding to the example porous media samples.

[0123] The predicted capillary pressure curves are generated using the alternative and enhanced pore-network modeling simulations, while the experimental data were measured using the mercury intrusion (MICP) and porous plate methods. All the values are adjusted to an interfacial tension of about 480 mN / m and a contact angle of about 140°. The pore space of the Bentheimersandstone is relatively homogeneous. As a result, the non-wetting phase can invade a large number of pores with a moderate increase in the capillary pressure, causing a quick reduction in the wetting phase saturation in the middle region of the capillary pressure curve. The predictions made by both PNM approaches reasonably follow the general trend of the measured capillary pressure data, indicating the extracted pore networks successfully capture the topology of the pore space. However, the alternative PNM approach tends to overestimate the capillary pressure curve for wetting-phase saturations above 0.1. This may be because the alternative pore-network models use simplified geometries, which limits the contribution of the corner space of a highly irregular pore element to the pore-level two-phase displacement process and tends to overestimate capillary entry pressures of pores. This problem is alleviated by applying the enhanced PNM platform, where the real geometries of the pore cross-sections are considered for the calculation. As a result, the predicted capillary pressure curve using the enhanced pore network shows a closer agreement with the experimental data in general.

[0124] Figure 15 is a graph of the predicted and measured capillary pressure curves for drainage in the Berea sandstone porous media sample. Figure 16 (d)- (f) is the pore fluid-occupancy maps predicted by the enhanced PNM approach for the Berea sample. The one or more processors perform primary drainage simulations in Berea sandstones and compare the predicted capillary pressure curves with the experimental data corresponding to the example porous media samples. The predicted capillary pressure curves are generated using the alternative and enhanced pore-network modeling simulations, while the experimental data were obtained from an image-based capillary pressure measurement using a micro-scale core-flooding apparatus and a porous membrane experiment. Contact angles of the oil-brine-Berea rock system used in the simulations are assigned a measured in-situ contact angle distribution. An interfacial tension value of about 41.3 mN / m is applied. In this example, the capillary pressure curve predicted by the enhanced PNM approach shows a good agreement with the experimental data, while the results generated from thealternative PNM consistently overestimate the measured values, as previously explained for the Bentheimer sample. Near the end of the oil-water primary drainage, where displacements primarily take place in small pores and crevices at high capillary pressures, any discrepancy between the steady-state miniature core-flooding experiment and the unsteady-state macro-scale porous membrane test could be amplified. Nevertheless, we note that simulation results from both the alternative and enhanced PNM approaches are reasonably aligned with the experimental data. Furthermore, the alternative PNM always predicts higher capillary pressure values compared to the enhanced PNM, as was similarly found for the Bentheimer sample. Overall, both PNM approaches show promising capability in predicting multiphase flow properties of homogenous sandstones, and improved performances are obtained from the enhanced PNM that better accounts for the effect of the pore geometries.

[0125] Figure 17 is a graph of the predicted capillary pressure curves for MICP test in the Fond du Lac dolostone (FLD) rock sample using the traditional and enhanced pore-network modeling approaches. Figure 18 is a graph of the predicted capillary pressure curves for MICP test in the Desert Sunset limestone rock (DSL) sample using the traditional and enhanced pore-network modeling approaches. After the validation of the enhanced PNM platform, we conduct MICP simulations in the conventional and enhanced pore networks of the FLD and DSL rock samples with more complicated pore structures. An interfacial tension of 480 mN / m and a contact angle of 140° are applied in both simulations.

[0126] Figure 19(a) is a visualization of the 3D invading-phase clusters of FLD at a saturation of 0.19. Figure 19(b) is a visualization of the 3D invading-phase clusters of FLD at a saturation of 0.41. Figure 19(c) is a visualization of the 3D invading-phase clusters of FLD at a saturation of 0.61. Figure 19(d) is a visualization of the 2D flu id -occupancy maps at a saturation of 0.19. Figure 19(e) is a visualization of the 2D fluid-occupancy maps of FLD at a saturation of 0.41. Figure 19(f) is a visualization of the 2D fluid-occupancy maps of FLD at a saturation of 0.61. The visualizations of the simulated non-wetting-phase clusters and fluid-occupancy maps in the real pore space during the primary drainage displacements in the FLD use the enhanced PNM approach.

[0127] Figure 20(a) is a visualization of the 3D invading-phase clusters of DSL at a saturation of 0.21. Figure 20(b) is a visualization of the 3D invading-phase clusters of DSL at a saturation of 0.42. Figure 20(c) is a visualization of the 3D invading-phase clusters of DSL at a saturation of 0.61. Figure 20(d) is a visualization of the 2D fluid-occupancy maps of DSL at a saturation of 0.21. Figure 20(e) is a visualization of the 2D fluid-occupancy maps of DSL at a saturation of 0.42. Figure 20(f) is a visualization of the 2D fluid-occupancy maps of DSL at a saturation of 0.61. The visualizations of the simulated non-wetting-phase clusters and fluid-occupancy maps during primary drainage in the DSL using the enhanced PNM approach.

[0128] The example capillary pressure curves are less smooth for the carbonates than sandstones. This may be due to the pore spaces of the two carbonates are more heterogeneous. As the non-wetting fluid is introduced into the example porous media sample, the invasion front can be hindered by small throat restrictions that connect some large void spaces of the main flow paths. Once the pressure of the invading fluid overcomes the capillary entry pressure of such throat restrictions, the invading fluid may quickly occupy a number of connected large pores, leading to a sudden increase in the capillary pressure curve. The irreducible wetting-phase saturation of the DSL is high due to a large number of isolated pore elements in the extracted pore network. Such isolated pore elements in the pore network may be connected through regions of unresolved micro-pores, whose sizes are below the resolution of the micro-CT images.

[0129] In some cases, in comparison with the homogeneous sandstone rock sample examples, the carbonate examples show more pronounced differences between the predicted capillary pressures obtained from the alternative and enhanced PNM approaches. This may be due, in part, to: 1) the shapes of the individual pores in carbonates are more complex; or 2) the spatial topology of the pore space in carbonates is more heterogeneous. As a result, the simplifiedrepresentation of the pore space with regular convex triangles or squares may lead to larger deviations in estimating the capillary entry pressure of the pores. In this case, the deviations in the computed pore-scale capillary entry pressures are more prone to generate distinct core-scale invasion patterns and fluid-occupancy maps, which ultimately leads to different capillary pressure curves. This corroborates the necessity of using the enhanced PNM with real pore geometries to predict the flow properties of complex rocks, such as carbonates.

[0130] Figure 21 is a graph of the Pearson correlation between the geometric features and the capillary entry pressure of pores. Figure 22 is a scatterplot matrix of the geometric features against the capillary entry pressure. The pores are extracted from the Bentheimer sandstone, Berea sandstone, Fond du Lac dolostone, and Desert Sunset limestone. The Pearson correlation coefficient is computed between each geometric feature and the capillary entry pressures predicted using our semi-analytical model.

[0131] It may be inferred that the order of the geometric features from having the strongest to the weakest impact on the capillary entry pressure is: hydraulic radius > mean distance > inscribed radius > equivalent diameter > circumradius > perimeter > sphericity > area > aspect ratio > roundness > convexity > shape factor > extent > solidity, though other orders are contemplated. In this example, geometric features having strong associations with the capillary entry pressure, the hydraulic radius, mean distance, inscribed radius, equivalent diameter, circumradius, and perimeter show negative correlations with the capillary entry pressure. This may be expected for a standard porous media sample, as the capillary entry pressure is usually smallerfor larger pore elements. In this example, several geometric features may demonstrate very strong correlations with each other. For example, the hydraulic radius has correlation coefficients of about 0.98, about 0.98, and about 0.97 against the inscribed radius, mean distance, and equivalent diameter, respectively. Additionally, the circumradius has correlation coefficients of about 0.98 against the perimeter and about 0.97 against the equivalent diameter. In this example, it may be reasonable to choose onerepresentative geometric feature from these similar features to quantify their impacts on the capillary entry pressure of pores.

[0132] To construct an inclusive ANN 1000 that can comprehensively predict capillary entry pressures of pores in different types of porous media, the one or more processors may randomly sample a total of between about 50,000 pores or more to about 250,000 pores or less (e.g., 100,000 pores), though other values are contemplated. In example, one or more processors may construct throat cross- sections from the four enhanced pore networks of the Bentheimer sandstone, Berea sandstone, FLD, and DSL. Their detailed geometric features may be analyzed according to a set of predefined feature definitions. Certain features may be discarded, and certain features given additional weight. In view of certain physical relationships in the example, the processors may replace the hydraulic radius and circumradius with their inverse as the input features to better inform the ANN 1000. Furthermore, as the target variable of the machine learning model, the capillary entry pressures of all the sampled pore elements may be computed using the semi-analytical method 500 following the fluid properties used in the previous MICP simulations. The final data set is split into training, validation, and test sets with a corresponding percentage of about 56.25%, about 18.75%, and about 25%.

[0133] After a randomized search on the hyperparameters space, the one or more processors may find that the ANN 1000 configuration consisting of about 2 hidden layers with about 7 neurons per hidden layer and having a learning rate of about 5.93 x io-4gives the optimal performance on capillary entry pressure predictions for this example. The mean training loss measured over each epoch and the mean validation loss obtained at the end of each epoch during the training of the ANN 1000 with the optimal hyperparameters. The training of the ANN uses an ensemble of pore and throat elements sampled from the Bentheimer sandstone, Berea sandstone, Fond du Lac dolostone, and Desert Sunset limestone.

[0134] Figure 23 is a graph between the predicted capillary entry pressures using the ANN 1000 and the corresponding capillary entry pressures. It may be observed that the ANN 1000 converges rapidly with small training and validationlosses after a few epochs, indicating no underfitting and overfitting in the constructed neural network. An ensemble of pore and throat elements sampled from the Bentheimer sandstone, Berea sandstone, Fond du Lac dolostone, and Desert Sunset limestone is used for training the ANN 1000. The regression score, R2, with respect to the test data set may be found to be about 0.995. Accordingly, the developed ANN 1000 may predict the drainage capillary entry pressure of pores formed in various types of porous media samples (e.g., sandstone and carbonate rock samples) by utilizing key pore geometric features, including the hydraulic radius, circumradius, sphericity, and area (in the order of high-to-low influence).

[0135] Aspects of the present disclosure provide a workflow that couples an enhanced PNM approach with a thermodynamically-consistent, semi-analytical method 500 for pore capillary entry pressure to simulate two-phase drainage flow in porous media. Compared to existing pore-scale flow models, this enhanced PNM approach has the advantage of describing pore-level, two-phase displacement processes using real pore geometries while maintaining a high computational efficiency. Additionally, an energy-based criterion for eliminating invalid arc meniscus (AMs) and a parallel computing design are provided to enhance the stability and the computational performance of the semi-analytical model. The improved workflow may be illustrated according to piston-like drainage examples. The predicted drainage capillary pressure curves from the example implementations of the work illustrate satisfactory improvements over the results generated from techniques known in the current art. Thus, aspects of the present disclosure are critical to faithfully predict multiphase flow and transport properties of the heterogeneous porous media, and to utilize a computationally efficient pore- scale model that considers both the intricate topology of the pore space and the detailed pore geometries. Accordingly, the techniques described herein may reduce porous media sample characterization errors to the benefit of all users seeking a more comprehensive understanding of any given porous media.Example Device

[0136] Figure 24 depicts aspects of an example porous media device 2400. In some aspect, the device 2400 includes one or more CPUs, one or more GPUs, or both as described above with respect to Figure 4.

[0137] The device 2400 includes a CPU processing system 2404 coupled to an image interface 2402 (e.g., a user interface or and / or an image generator such as a commercial micro-CT scanner). The CPU processing system 2404 may be configured to perform processing functions for the device 2400, including multiscale imaging of porous media generated by the device 2400.

[0138] The CPU processing system 2404 includes one or more processors 2410. The one or more processors 2410 are coupled to a computer-readable medium / memory 2412 via a bus. The one or more processors 2410 and the computer-readable medium / memory 2412 may communicate with the one or more processor 2414 and the computer-readable medium / memory 2416 of the GPU processing system 2406 via a message passing interface (MPI) 2408. In certain aspects, the computer-readable medium / memory 2412 is configured to store instructions (e.g., computer-executable code) that when executed by the one or more processors 2410, cause the one or more processors 2410 to perform the method 400 described with respect to Figure 4, or any aspect related to it. Note that reference to a processor performing a function of device 2400 may include one or more processors performing that function of device 2400.

[0139] In the depicted example, computer-readable medium / memory 2412 stores code 2430, 2432, 2434, 2436, 2438, and 2440 (e.g., executable instructions) for performing certain aspects of the present disclosure. Processing of the code 2430-2440 may cause the imaging device 2400 to perform the method 400 described with respect to Figure 4, or any aspect related to it.

[0140] The one or more processors 2410 include circuitry configured to implement (e.g., execute) the code stored in the computer-readable medium / memory 2412, including circuitry 2418, 2420, 2422, 2424, 2426 and 2428 for performing certain aspects of the present disclosure. Processing with circuitry 2418-2428 may cause the imaging device 2400 to perform the method 400 described with respect to Figure 4, or any aspect related to it.

[0141] The device 2400 includes a GPU processing system 2406. The GPU processing system 2406 may be configured to perform processing functions for the device 2400, including multiscale imaging of porous media generated by the device 2400.

[0142] The GPU processing system 2406 includes one or more processors 2414. The one or more processors 2414 are coupled to a computer-readable medium / memory 2416 via a bus. The one or more processors 2414 and the computer-readable medium / memory 2416 may communicate with the one or more processor 2410 and the computer-readable medium / memory 2412 of the CPU processing system 2404 via an MPI 2408. In certain aspects, the computer- readable medium / memory 2416 is configured to store instructions (e.g., computer- executable code) that when executed by the one or more processors 2414, cause the one or more processors 2414 to perform the method 400 described with respect to Figure 4, or any aspect related to it. Note that reference to a processor performing a function of device 2400 may include one or more processors performing that function of device 2400.

[0143] In the depicted example, computer-readable medium / memory 2416 stores code 2452, 2454, 2456, 2458, and 2460 (e.g., executable instructions) for performing certain aspects of the present disclosure. Processing of the code 2452- 2460 may cause the imaging device 2400 to perform the method 400 described with respect to Figure 4, or any aspect related to it.

[0144] The one or more processors 2414 include circuitry configured to implement (e.g., execute) the code stored in the computer-readable medium / memory 2416, including circuitry 2442, 2444, 2446, 2448, and 2450 for performing certain aspects of the present disclosure. Processing with circuitry 2442-2450 may cause the imaging device 2400 to perform the method 400 described with respect to Figure 4, or any aspect related to it.Embodiments Listing

[0145] Implementation examples are described in the following numbered embodiments:

[0146] Clause 1 : A method for predicting the physical and chemical media properties of a porous media by one or more processing units, comprising: performing a pore-network modeling (PNM); applying a Mayer-Stowe-Princen (MS-P) semi-analytical theory of drainage to the PNM, comprising: constructing pore boundaries from an image of a pore cross-section; identifying geometrically and thermodynamically-favorable arc meniscus (AM) from aa plurality of AMs, the geometrically and thermodynamically-favorable AM having a contact angle and a curvature radius; determining an optimal fluid configuration from the geometrically and thermodynamically-favorable AMs.

[0147] Clause 2: The method of clause 1 , wherein constructing pore boundaries comprises: obtaining a high-resolution image of the pore cross section; extracting a pore boundary from the high-resolution image of the pore cross section utilizing the one or more processors having a computer vision and image analysis library to extract a boundary point of the high-resolution image of the pore cross section; smoothing the pore boundary using a moving-average filter, wherein the one or more processors apply the moving-average filter along the boundary point to smooth out artificial discontinuities of the pore boundary.

[0148] Clause 3: The method of clause 1 , wherein identifying the geometrically and thermodynamically favorable AMs comprises: removing AMs not forming a valid contact angle; removing AMs having a diverging curvatures; removing AMs outside of the pore boundary; and removing AMs that intersect with other AMs.

[0149] Clause 4: The method of clause 1 , wherein determining an optimal fluid configuration comprises: identifying phase regions;identifying an optimal fluid configuration using the one or more processors; updating an AM curvature; performing a query regarding whether the AM curvature has converged; and when the AM curvature has not converged, repeating applying a Mayer- Stowe-Princen (MS-P) semi-analytical theory of drainage to the PNM.

[0150] Clause 5: The method of clause 1 , further comprising: training an artificial neural network (ANN) utilizing the MS-P theory based semi-analytical theory of drainage; and estimating the physical and chemical media properties of pore using the ANN.

[0151] Clause 6: An apparatus, comprising: a memory comprising executable instructions; and a processor configured to execute the executable instructions and cause the apparatus to perform a method, the method comprising: performing a pore-network modeling (PNM); applying a Mayer-Stowe-Princen (MS-P) semi-analytical theory of drainage to the PNM, comprising: constructing pore boundaries from an image of a pore cross-section; identifying geometrically and thermodynamically-favorable arc meniscus (AM) from aa plurality of AMs, the geometrically and thermodynamically-favorable AM having a contact angle and a curvature radius; determining an optimal fluid configuration from the geometrically and thermodynamically-favorable AMs.

[0152] Clause 7: The apparatus of clause 6, wherein constructing pore boundaries comprises: obtaining a high-resolution image of the pore cross section; extracting a pore boundary from the high-resolution image of the pore cross section utilizing the one or more processors having a computer vision and image analysis library to extract a boundary point of the high-resolution image of the pore cross section;smoothing the pore boundary using a moving-average filter, wherein the one or more processors apply the moving-average filter along the boundary point to smooth out artificial discontinuities of the pore boundary.

[0153] Clause 8: The apparatus of clause 6, wherein identifying the geometrically and thermodynamically favorable AMs comprises: removing AMs not forming a valid contact angle; removing AMs having a diverging curvatures; removing AMs outside of the pore boundary; and removing AMs that intersect with other AMs.

[0154] Clause 9: The apparatus of clause 6, wherein determining an optimal fluid configuration comprises: identifying phase regions; identifying an optimal fluid configuration using the one or more processors; updating an AM curvature; performing a query regarding whether the AM curvature has converged; and when the AM curvature has not converged, repeating applying a Mayer- Stowe-Princen (MS-P) semi-analytical theory of drainage to the PNM.

[0155] Clause 10: The apparatus of clause 6, further comprising: training an artificial neural network (ANN) utilizing the MS-P theory based semi-analytical theory of drainage; and estimating the physical and chemical media properties of pore using the ANN.

[0156] Clause 11 : A non-transitory computer-readable medium comprising executable instructions that, when executed by a processor of an apparatus, cause the apparatus to perform a method, the method comprising: performing a pore-network modeling (PNM); applying a Mayer-Stowe-Princen (MS-P) semi-analytical theory of drainage to the PNM, comprising: constructing pore boundaries from an image of a pore cross-section;identifying geometrically and thermodynamically-favorable arc meniscus (AM) from a plurality of AMs, the geometrically and thermodynamically-favorable AM having a contact angle and a curvature radius; determining an optimal fluid configuration from the geometrically and thermodynamically-favorable AMs.

[0157] Clause 12: The non-transitory computer-readable medium of clause 11 , wherein constructing pore boundaries comprises: obtaining a high-resolution image of the pore cross section; extracting a pore boundary from the high-resolution image of the pore cross section utilizing the one or more processors having a computer vision and image analysis library to extract a boundary point of the high-resolution image of the pore cross section; smoothing the pore boundary using a moving-average filter, wherein the one or more processors apply the moving-average filter along the boundary point to smooth out artificial discontinuities of the pore boundary.

[0158] Clause 13: The non-transitory computer-readable medium of clause 11 , wherein identifying the geometrically and thermodynamically favorable AMs comprises: removing AMs not forming a valid contact angle; removing AMs having a diverging curvatures; removing AMs outside of the pore boundary; and removing AMs that intersect with other AMs.

[0159] Clause 14: The non-transitory computer-readable medium of clause 11 , wherein determining an optimal fluid configuration comprises: identifying phase regions; identifying an optimal fluid configuration using the one or more processors; updating an AM curvature; performing a query regarding whether the AM curvature has converged; and when the AM curvature has not converged, repeating applying a Mayer- Stowe-Princen (MS-P) semi-analytical theory of drainage to the PNM.

[0160] Clause 15: The non-transitory computer-readable medium of clause 11 , further comprising: training an artificial neural network (ANN) utilizing the MS-P theory based semi-analytical theory of drainage; and estimating the physical and chemical media properties of pore using the ANN.Additional Considerations

[0161] The preceding description is provided to enable any person skilled in the art to practice the various aspects described herein. The examples discussed herein are not limiting of the scope, applicability, or aspects set forth in the claims. Various modifications to these aspects will be readily apparent to those skilled in the art, and the general principles defined herein may be applied to other aspects. For example, changes may be made in the function and arrangement of elements discussed without departing from the scope of the disclosure. Various examples may omit, substitute, or add various procedures or components as appropriate. For instance, the methods described may be performed in an order different from that described, and various actions may be added, omitted, or combined. Also, features described with respect to some examples may be combined in some other examples. For example, an apparatus may be implemented or a method may be practiced using any number of the aspects set forth herein. In addition, the scope of the disclosure is intended to cover such an apparatus or method that is practiced using other structure, functionality, or structure and functionality in addition to, or other than, the various aspects of the disclosure set forth herein. It should be understood that any aspect of the disclosure disclosed herein may be embodied by one or more elements of a claim.

[0162] As used herein, a phrase referring to “at least one of’ a list of items refers to any combination of those items, including single members. As an example, “at least one of: a, b, or c” is intended to cover a, b, c, a-b, a-c, b-c, and a-b-c, as well as any combination with multiples of the same element (e.g., a-a, a-a-a, a-a-b, a- a-c, a-b-b, a-c-c, b-b, b-b-b, b-b-c, c-c, and c-c-c or any other ordering of a, b, and c). The singular forms “a,” “an,” and “the” include plural referents, unless thecontext clearly dictates otherwise. Within a claim, reference to an element in the singular is not intended to mean “one and only one” unless specifically so stated, but rather “one or more.” Unless specifically stated otherwise, the term “some” refers to one or more.

[0163] As used herein, the term “determining" encompasses a wide variety of actions. For example, "determining” may include calculating, computing, processing, deriving, investigating, updating, looking up (e.g., looking up in a table, a database or another data structure), ascertaining and the like. Also, “determining" may include receiving (e.g., receiving information), accessing (e.g., accessing data in a memory) and the like. Also, "determining” may include resolving, selecting, simulating, choosing, establishing, and the like.

[0164] The methods disclosed herein include one or more operations or actions for achieving the methods. The method operations and / or actions may be interchanged with one another without departing from the scope of the claims. In other words, unless a specific order of operations or actions is specified, the order and / or use of specific operations and / or actions may be modified without departing from the scope of the claims. Further, the various operations of methods described above may be performed by any suitable means capable of performing the corresponding functions. The means may include various hardware and / or software component(s) and / or module(s), including, but not limited to a circuit, an application specific integrated circuit (ASIC), or processor. Generally, where there are operations illustrated in the figures, those operations may have corresponding counterpart means-plus-function components with similar numbering.

[0165] When the word “approximately" or “about” are used, this term may mean that there may be a variance in value of up to ±10%, of up to 5%, of up to 2%, of up to 1%, of up to 0.5%, of up to 0.1%, or up to 0.01 %.

[0166] Ranges may be expressed as from about one particular value to about another particular value, inclusive. When such a range is expressed, it is to be understood that another embodiment is from the one particular value to the other particular value, along with all particular values and combinations thereof within the range.

[0167] As used, terms such as “first” and “second” are arbitrarily assigned and are merely intended to differentiate between two or more components of a system, an apparatus, or a composition. It is to be understood that the words “first” and “second” serve no other purpose and are not part of the name or description of the component, nor do they necessarily define a relative location or position of the component. Furthermore, it is to be understood that that the mere use of the term “first” and “second” does not require that there be any “third” component, although that possibility is envisioned under the scope of the various embodiments described.

[0168] The following claims are not intended to be limited to the aspects shown herein, but are to be accorded the full scope consistent with the language of the claims. Within a claim, reference to an element in the singular is not intended to mean “one and only one” unless specifically so stated, but rather "one or more.” Unless specifically stated otherwise, the term “some” refers to one or more. No claim element is to be construed under the provisions of 35 U.S.C. §112(f) unless the element is expressly recited using the phrase “means for” or, in the case of a method claim, the element is recited using the phrase “step for.” All structural and functional equivalents to the elements of the various aspects described throughout this disclosure that are known or later come to be known to those of ordinary skill in the art are expressly incorporated herein by reference and are intended to be encompassed by the claims. Moreover, nothing disclosed herein is intended to be dedicated to the public regardless of whether such disclosure is explicitly recited in the claims.

[0169] Unless defined otherwise, all technical and scientific terms used have the same meaning as commonly understood by one of ordinary skill in the art to which these systems, apparatuses, methods, processes and compositions belong.

[0170] The following claims are not intended to be limited to the embodiments provided but rather are to be accorded the full scope consistent with the language of the claims.

Claims

AMENDED CLAIMS received by the International Bureau on 6 January 2025 (06.01.2025)1. A method for predicting physical and chemical media properties of a porous media by one or more processors, comprising: imaging a porous media using a micro-CT scanner to form a 3D image of the porous media; performing a pore-network extraction (PNE) to extract a 2D image of a pore cross- section of a pore within the porous media ; applying a Mayer-Stowe-Princen (MS-P) semi-analytical theory of drainage to the 2D image of the pore cross-section , comprising: constructing pore boundaries from the 2D image of the pore cross-section; identifying geometrically and thermodynamically-favorable arc meniscus (AM) from a plurality of AMs, the geometrically and thermodynamically-favorable AM having a contact angle and a curvature radius; and determining aa fluid configuration from the geometrically and thermodynamically-favorable AMs; coupling a pore-network modeling (PNM) with the MS-P semi-analytical theory of drainage; and performing the PNM enhanced with the MS-P semi-analytical theory of drainage to form a 3D pore network model using the fluid configuration.

2. The method of claim 1, wherein constructing pore boundaries comprises: obtaining a high-resolution image of the pore cross section; extracting a pore boundary from the high-resolution image of the pore cross section utilizing the one or more processors having a computer vision and image analysis library to extract a boundary point of the high-resolution image of the pore cross section; and smoothing the pore boundary using a moving-average filter, wherein the one or more processors apply the moving-average filter along the boundary point to smooth out artificial discontinuities of the pore boundary.

3. The method of claim 1, wherein identifying the geometrically and thermodynamically favorable AMs comprises: removing AMs not forming a valid contact angle;removing AMs having a diverging curvatures; removing AMs outside of the pore boundary; and removing AMs that intersect with other AMs.

4. The method of claim 1, wherein determining an optimal fluid configuration comprises: identifying phase regions; identifying an optimal fluid configuration using the one or more processors; updating an AM curvature; performing a query regarding whether the AM curvature has converged; and when the AM curvature has not converged, repeating applying the MS-P semi- analytical theory of drainage to the PNM.

5. The method of claim 1, further comprising: training an artificial neural network (ANN) utilizing the MS-P semi-analytical theory of drainage; and estimating the physical and chemical media properties of the pore using the ANN.

6. An apparatus, comprising: a memory comprising executable instructions; and a processor configured to execute the executable instructions and cause the apparatus to perform a method, the method comprising: imaging a porous media using a micro-CT scanner to form a 3D image of the porous media; performing a pore-network extraction (PNE) to extract a 2D image of a pore cross- section of a pore within the porous media ; applying a Mayer-Stowe-Princen (MS-P) semi-analytical theory of drainage to the 2D image of the pore cross-section , comprising: constructing pore boundaries from the 2D image of the pore cross-section; identifying geometrically and thermodynamically-favorable arc meniscus (AM) from a plurality of AMs, the geometrically and thermodynamically-favorable AM having a contact angle and a curvature radius; and determining aa fluid configuration from the geometrically and thermodynamically-favorable AMs;coupling a pore-network modeling (PNM) with the MS-P semi-analytical theory of drainage; and performing the PNM enhanced with the MS-P semi-analytical theory of drainage to form a 3D pore network model using the fluid configuration.

7. The apparatus of claim 6, wherein constructing pore boundaries comprises: obtaining a high-resolution image of the pore cross section; extracting a pore boundary from the high-resolution image of the pore cross section utilizing the one or more processors having a computer vision and image analysis library to extract a boundary point of the high-resolution image of the pore cross section; and smoothing the pore boundary using a moving-average filter, wherein the one or more processors apply the moving-average filter along the boundary point to smooth out artificial discontinuities of the pore boundary.

8. The apparatus of claim 6, wherein identifying the geometrically and thermodynamically favorable AMs comprises: removing AMs not forming a valid contact angle; removing AMs having a diverging curvatures; removing AMs outside of the pore boundary; and removing AMs that intersect with other AMs.

9. The apparatus of claim 6, wherein determining an optimal fluid configuration comprises: identifying phase regions; identifying an optimal fluid configuration using the one or more processors; updating an AM curvature; performing a query regarding whether the AM curvature has converged; and when the AM curvature has not converged, repeating applying the MS-P semi- analytical theory of drainage to the PNM.

10. The apparatus of claim 6, further comprising: training an artificial neural network (ANN) utilizing the MS-P semi-analytical theory of drainage; and estimating the physical and chemical media properties of the pore using the ANN.

11. A non-transitory computer-readable medium comprising executable instructions that, when executed by a processor of an apparatus, cause the apparatus to perform a method, the method comprising: imaging a porous media using a micro-CT scanner to form a 3D image of the porous media; performing a pore-network extraction (PNE) to extract a 2D image of a pore cross- section of a pore within the porous media ; applying a Mayer-Stowe-Princen (MS-P) semi-analytical theory of drainage to the 2D image of the pore cross-section , comprising: constructing pore boundaries from the 2D image of the pore cross-section; identifying geometrically and thermodynamically-favorable arc meniscus (AM) from a plurality of AMs, the geometrically and thermodynamically-favorable AM having a contact angle and a curvature radius; and determining aa fluid configuration from the geometrically and thermodynamically-favorable AMs; coupling a pore-network modeling (PNM) with the MS-P semi-analytical theory of drainage; and performing the PNM enhanced with the MS-P semi-analytical theory of drainage to form a 3D pore network model using the fluid configuration.

12. The non-transitory computer-readable medium of claim 11, wherein constructing pore boundaries comprises: obtaining a high-resolution image of the pore cross section; extracting a pore boundary from the high-resolution image of the pore cross section utilizing the one or more processors having a computer vision and image analysis library to extract a boundary point of the high-resolution image of the pore cross section; and smoothing the pore boundary using a moving-average filter, wherein the one or more processors apply the moving-average filter along the boundary point to smooth out artificial discontinuities of the pore boundary.

13. The non-transitory computer-readable medium of claim 11, wherein identifying the geometrically and thermodynamically favorable AMs comprises: removing AMs not forming a valid contact angle;removing AMs having a diverging curvatures; removing AMs outside of the pore boundary; and removing AMs that intersect with other AMs.

14. The non-transitory computer-readable medium of claim 11, wherein determining an optimal fluid configuration comprises: identifying phase regions; identifying an optimal fluid configuration using the one or more processors; updating an AM curvature; performing a query regarding whether the AM curvature has converged; and when the AM curvature has not converged, repeating applying the MS-P semi- analytical theory of drainage to the PNM.

15. The non-transitory computer-readable medium of claim 11, further comprising: training an artificial neural network (ANN) utilizing the MS-P semi-analytical theory of drainage; and estimating the physical and chemical media properties of the pore using the ANN.

Images