Methods and systems for simulating multiphase flow through porous media

Abstract

The present disclosure provides for generating multiphase flow properties of porous media based on one or more input parameters. For instance, the multiphase flow properties may be a capillary pressure-saturation relationship for the porous media. The one or more input parameters include an interfacial tension along an interface between a wetting fluid and a non-wetting fluid, a contact angle between the interface and a pore wall of the porous media, and a pore throat size. The pore throat size is based on subparameters including a saturation of the wetting fluid, a saturation of the non-wetting fluid, a porosity of the porous media, and an orientation angle between a representative pore body size and a representative pore throat size.

Claims

1. A non-transitory, computer-readable medium storing instructions which, when performed by a processor, cause the processor to: responsive to receiving one or more input parameters, generate a pressure-saturation relationship for a porous media thereby simulating capillary pressure-saturation behavior from one or more sample sizes based on the generated pressure-saturation relationship to reduce computational cost, wherein the one or more input parameters include an interfacial tension along an interface between a wetting fluid and a non-wetting fluid, a contact angle between the interface and a pore wall of the porous media, and a pore throat size, wherein the pore throat size is based on subparameters including a saturation of the wetting fluid, a saturation of the non-wetting fluid, a porosity of the porous media, and an orientation angle between a representative pore body size and a representative pore throat size, and wherein the representative pore body size (r.sub.b), the representative pore throat size (r.sub.t), and the orientation angle (β) are related according to the below relationship, $r_b=\frac{r_t}{1-\tan \left(\beta \right)}.$

2. The non-transitory, computer-readable medium of claim 1, wherein the pressure-saturation relationship is generated according to the below relationship of input parameters, $p_c\left(s_w\right)=\frac{2\sigma \cos \left(\theta \right)}{{\hat{R}}_t}$ wherein p.sub.c is capillary pressure, s.sub.w is the saturation of the wetting fluid, σ is the interfacial tension, θ is the contact angle, and {circumflex over (R)}.sub.t is the pore throat size.

3. The non-transitory, computer-readable medium of claim 1, wherein the pore throat size {circumflex over (R)}.sub.t is based on the below relationship of subparameters, ${\hat{R}}_t=\frac{{\stackrel{\_}{R}}_t\left(s_n\right)}{\left(S_w^a+{\left(1-s_w\right)}^2\right)\exp \left(s_w\right)},\mathrm{wherein}$ ${\stackrel{\_}{R}}_t\left(s_n\right)={\stackrel{\_}{R}}_t\left(s_{n-1}\right)-\frac{\sqrt[3]{\varnothing {s_{n-1}\left(1-\tan \beta \right)}^3{r}_g^3}}{3}\times \hspace{1em}\left[\frac{\frac{1}{2}+\frac{3}{2}{\left(1-\tan \beta \right)}^3}{{\left[s_{n-1}\left[\frac{1}{2}+\frac{3}{2}{\left(1-\tan \beta \right)}^3\right]-\varnothing s_{n-1}\left[\frac{1}{2}+\frac{3}{2}{\left(1-\tan \beta \right)}^3\right]\right]}^{\frac{4}{3}}}\right]\times \left[s_n-s_{n-1}\right],$ wherein s.sub.w is the saturation of the wetting fluid, s.sub.n is the saturation of the non-wetting fluid, β is the orientation angle, Ø is the porosity, r.sub.g is a representative grain size, and a is the rugosity.

4. The non-transitory, computer-readable medium of claim 1, wherein a pore body size {circumflex over (R)}.sub.b is related to the pore throat size {circumflex over (R)}.sub.t according to the below relationship, ${\hat{R}}_b=\frac{{\hat{R}}_t}{1-\tan \left(\beta \right)}$ wherein β is the orientation angle.

5. The non-transitory, computer-readable medium of claim 1, wherein the pressure-saturation relationship is generated according to the below relationship of input parameters, $p_c\left(s_w\right)=\frac{2\sigma \cos \left(\theta \right)}{{\hat{R}}_b}$ wherein p.sub.c is capillary pressure, s.sub.w is the saturation of the wetting fluid, σ is the interfacial tension, and θ is the contact angle.

6. The non-transitory, computer-readable medium of claim 1, wherein the orientation angle is based on a mean pore body size and a mean pore throat size.

7. The non-transitory, computer-readable medium of claim 1, wherein one of the orientation angle and the porosity is extracted from one or more images.

8. The non-transitory, computer-readable medium of claim 1, wherein the subparameters further include a representative grain size and a rugosity of the porous media.

9. A system for forecasting multiphase flow properties through porous media, the system comprising: a processor; and a memory storing instructions which, when executed by the processor, cause the processor to: responsive to receiving one or more input parameters, generate a pressure-saturation relationship for a porous media thereby simulating capillary pressure-saturation behavior from one or more sample sizes based on the generated pressure-saturation relationship to reduce computational cost, wherein the one or more input parameters include an interfacial tension along an interface between a wetting fluid and a non-wetting fluid, a contact angle between the interface and a pore wall of the porous media, and a pore throat size, and wherein a pore body size {circumflex over (R)}.sub.b is related to the pore throat size {circumflex over (R)}.sub.t according to the below relationship, ${\hat{R}}_b=\frac{{\hat{R}}_t}{1-\tan \left(\beta \right)}$ wherein β is an orientation angle.

10. The system of claim 9, wherein the pore throat size is based on subparameters including a saturation of the wetting fluid, a saturation of the non-wetting fluid, a porosity of the porous media, the orientation angle between a representative pore body size and a representative pore throat size, a representative grain size, and a rugosity of the porous media.

11. The system of claim 9, wherein one of the porosity and the orientation angle is determined through one or more images.

12. The system of claim 9, wherein the pressure-saturation relationship is generated according to the below relationship of input parameters, $p_c\left(s_w\right)=\frac{2\sigma \cos \left(\theta \right)}{{\hat{R}}_t}$ wherein p.sub.c is capillary pressure, s.sub.w is the saturation of the wetting fluid, σ is the interfacial tension, θ is the contact angle, and {circumflex over (R)}.sub.t is the pore throat size.

13. The system of claim 10, wherein the pore throat size {circumflex over (R)}.sub.t is based on the below relationship of subparameters, ${\hat{R}}_t=\frac{{\stackrel{\_}{R}}_t\left(s_n\right)}{\left(S_w^a+{\left(1-s_w\right)}^2\right)\exp \left(s_w\right)},\mathrm{wherein}$ ${\stackrel{\_}{R}}_t\left(s_n\right)={\stackrel{\_}{R}}_t\left(s_{n-1}\right)-\frac{\sqrt[3]{\varnothing {s_{n-1}\left(1-\tan \beta \right)}^3{r}_g^3}}{3}\times \hspace{1em}\left[\frac{\frac{1}{2}+\frac{3}{2}{\left(1-\tan \beta \right)}^3}{{\left[s_{n-1}\left[\frac{1}{2}+\frac{3}{2}{\left(1-\tan \beta \right)}^3\right]-\varnothing s_{n-1}\left[\frac{1}{2}+\frac{3}{2}{\left(1-\tan \beta \right)}^3\right]\right]}^{\frac{4}{3}}}\right]\times \left[s_n-s_{n-1}\right],$ wherein s.sub.w is the saturation of the wetting fluid, s.sub.n is the saturation of the non-wetting fluid, β is the orientation angle, Ø is the porosity, r.sub.g is the representative grain size, and a is the rugosity.

14. The system of claim 9, wherein the pressure-saturation relationship is generated according to the below relationship of input parameters, $p_c\left(s_w\right)=\frac{2\sigma \cos \left(\theta \right)}{{\hat{R}}_b}$ wherein p.sub.c is capillary pressure, s.sub.w is the saturation of the wetting fluid, σ is the interfacial tension, and θ is the contact angle.

15. The system of claim 9, wherein the orientation angle is based on a mean pore body size and a mean pore throat size.

16. A method for forecasting multiphase flow properties through porous media, the method comprising: determine one or more input parameters, wherein the one or more input parameters include a porosity of a porous media, an orientation angle between a representative pore body size and a representative pore throat size, and a representative grain size; and input the one or more parameters to generate a pressure-saturation relationship for the porous media thereby simulating capillary pressure-saturation behavior from one or more sample sizes based on the generated pressure-saturation relationship to reduce computational cost, wherein the pressure-saturation relationship is based on an interfacial tension along an interface between the wetting fluid and the non-wetting fluid, a contact angle between the interface and a pore wall of the porous media, the porosity, the orientation angle, and the representative grain size, wherein a pore body size {circumflex over (R)}.sub.b is related to the pore throat size {circumflex over (R)}.sub.t according to the below relationship, ${\hat{R}}_b=\frac{{\hat{R}}_t}{1-\tan \left(\beta \right)}$ wherein β is the orientation angle.

17. The method of claim 16, wherein at least one of the porosity or the orientation angle is determined through one or more images.

18. The method of claim 16, wherein the method further comprises determining the representative grain size such that the pore throat size decreases or is constant as a saturation of the wetting fluid increases.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1 illustrates a graph showing capillary pressure-saturation trends of porous media.

(2) FIG. 2 illustrates a representative elementary volume (REV) of a porous media represented as a hypothesized homogenous porous medium.

(3) FIG. 3A illustrates a hypothesized homogenous REV.

(4) FIG. 3B illustrates converging-diverging channels representing a pore geometry.

(5) FIG. 4 illustrates a graph showing pore throat size progressively decreasing with the saturation of a wetting fluid.

(6) FIG. 5 illustrates an example schematic of a roughness coefficient as a function of the saturation of a wetting fluid under drainage conditions.

(7) FIG. 6 illustrates example trends of a roughness coefficient as a function of the saturation of a wetting fluid in response to variations in a rugosity parameter.

(8) FIG. 7 illustrates an example schematic showing variation in a roughness coefficient as a result of the saturation of a wetting fluid for the corner flow mechanism.

(9) FIG. 8 illustrates an example method for forecasting multiphase flow properties through porous media under drainage conditions, according to an aspect of the present disclosure.

(10) FIG. 9 illustrates an example method for forecasting multiphase flow properties through porous media under imbibition conditions, according to an aspect of the present disclosure.

(11) FIG. 10 illustrates an example system for forecasting multiphase flow properties through porous media, according to an aspect of the present disclosure.

(12) FIG. 11 illustrates the results of a sensitivity analysis performed for capillary pressure-saturation relationships generated by the provided methods.

(13) FIG. 12 illustrates a comparison between the provided method and conventional methods for forecasting multiphase flow properties of various Berea sandstone samples with differing grain textural properties under drainage conditions.

(14) FIG. 13 illustrates a comparison between the provided method and conventional methods for forecasting multiphase flow properties of various Berea sandstone samples with differing grain textural properties under imbibition conditions.

(15) FIG. 14 illustrates a comparison between the provided method and data involving displacement of water by air in sand packs.

(16) FIG. 15 illustrates a comparison between the provided method and pore network modeling.

(17) FIG. 16 illustrates a comparison between the provided method and pore network modeling.

DETAILED DESCRIPTION

(18) Capillary pressure and relative permeability measurements are an integral part of special core analysis (SCAL), a laboratory procedure for conducting flow experiments on core plugs taken from a petroleum reservoir, which the oil and gas industry greatly relies on. Capillary pressure measurements are required to determine the thickness of the water-oil transition zone, and to perform some displacement calculations. Reservoir engineers implement the capillary pressure measurements in simulators to determine the amount of hydrocarbons as well as the flowing capacity of fluids in a petroleum reservoir. Relative permeability measurements are required to predict flow in the reservoir and to model reservoir displacement processes. Despite the importance of capillary pressure and relative permeability measurements, conventional laboratory and pore-scale simulation techniques used to measure capillary pressure curves of core samples are expensive, tedious, time-consuming and prone to error.

(19) Accordingly, the presently disclosed methods and system provide for a new analytical methodology that can provide a reliable forecast of capillary pressure-saturation relationships and relative permeability of porous media under static and dynamic conditions, recovery efficiency of fluids, and other properties of multiphase flow through porous media. The provided analytical methodology may represent pore space with a simple idealized geometry. The provided analytical methodology also requires very few input parameters to generate a capillary pressure-saturation relationship for a porous media that, in some instances, can be determined from pore-scale images of reservoir rocks. In other instances, the provided analytical methodology may utilize average pore properties rather than complete pore size distribution, and thus does not rely on input parameters extracted from images.

(20) As such, a capillary pressure-saturation relationship for a porous media may be generated using the provided analytical methodology by solving simple analytical calculations with a data analysis tool (e.g., Microsoft® Excel). For instance, input parameters may be introduced into the equations in the method described below, and a data analysis tool may solve the equations consistent with the provided method to generate the capillary pressure-saturation relationship. In some aspects of the present disclosure, the provided analytical methodology may be implemented into a software application configured to generate capillary pressure-saturation relationships for porous media upon the input of the described parameters.

(21) Thus, the presently disclosed method's analytical formulation allows it to simulate the capillary pressure-saturation relationship of porous media faster than conventional methods and may be free from numerical instabilities associated with numerical solvers. For instance, the accuracy of the provided systems and methods may depend solely on the input data. Moreover, the analytical nature makes the presently disclosed method simple to implement, highly reproducible, and flexible enough to be capable of predicting capillary pressure-saturation behavior from any sample size (e.g., pore- to core-scale). Additionally, the provided method is computationally inexpensive to solve, and thus may have an advantage of being amendable towards integration within existing continuum-scale modelling frameworks.

(22) Therefore, the provided analytical methodology may provide the oil and gas industry an inexpensive, fast, and accurate estimation of capillary pressure-saturation data. The provided methodology may also reduce the number of required laboratory experiments and may facilitate the estimation of capillary pressure-saturation data from un-cored sections of a petroleum reservoir (e.g., using drill cuttings). In addition to the oil and gas industry, the provided analytical methodology may provide a valuable tool towards the study of a broad range of porous media in a variety of other industries (e.g., energy and environment, textile, agriculture, etc.) as the methodology can reduce the need for expensive experiments and time-consuming numerical simulations to determine the multiphase flow properties (e.g., capillary pressure-saturation characteristics) of a given sample.

(23) Throughout this disclosure, reference is made to a number of equations and symbols. Unless stated otherwise, each of the symbols in the various equations has the definition given in Table 1 below. The subscripts n and w refer to non-wetting fluid and wetting fluid, respectively.

(24) TABLE-US-00001 TABLE 1 Nomenclature Symbol Description i i = n (non-wetting fluid) and i = w (wetting fluid) p.sub.c Capillary Pressure p.sub.n Pressure of the non-wetting fluid p.sub.w Pressure of the wetting fluid Ø.sub.REV Porosity of hypothesized REV Ø Porosity of porous media under investigation V.sub.l* Fluid volumes within the REV V* Total volume of REV s.sub.t Saturation of fluid in porous media under investigation r.sub.g Grain size of hypothesized REV r.sub.b Pore Body size of hypothesized REV r.sub.t Pore throat size of hypothesized REV N.sub.g Number of pore bodies in hypothesized REV N.sub.t Number of pore throats in hypothesized REV N.sub.G Number of grains in hypothesized REV s.sub.i* Saturation of fluids in hypothesized REV x.sub.b Correlation factor between N.sub.g and N.sub.G x.sub.t Correlation factor between N.sub.t and N.sub.G β Orientation angle between r.sub.t and r.sub.b of hypothesized REV R.sub.i Pore throat size of porous media under investigation s.sub.wr Irreducible wetting phase saturation in porous media under investigation s.sub.n−1 Previous saturation of non-wetting fluid in porous media under investigation s.sub.w−1 Previous saturation of wetting fluid in porous media under investigation t Number of capillary pressure-saturation data points r.sub.min Minimum pore throat size of porous media under investigation f.sub.r(s.sub.w) Roughness coefficient {hacek over (R)}.sub.t The effective pore throat size of porous media under investigation a The rugosity parameter for computation of f.sub.r(s.sub.w) R.sub.b The effective pore body size of porous media under investigation P.sub.d The displacement pressure V.sub.b Fractional bulk mercury saturation V.sub.∞ Fractional bulk mercury saturation at infinite p.sub.c F.sub.g The pore geometric factor

(25) One concept important to modelling multiphase flow in porous media is Representative Elementary Volume (“REV”), which is the smallest volume over which a measurement can be made that will yield a value representative of the whole. According to REV concepts, there exists a scale at which heterogeneities in measured properties within a porous media become statistically homogenous and immune to boundary effects to form an effective continuum medium. Hill, R., Elastic properties of reinforced solids: some theoretical principles. J. Mech. Phys. Solids 11 (5), 357-372 (1963). The minimum volume at which this statistical homogenization occurs at is the minimum representative element (e.g., REV) of the porous medium for a given property. As shown in the capillary pressure-saturation trends of FIG. 1, REV exists for the capillary pressure-saturation relationship in porous media. Joekar-Niasar, V., et al., Non-equilibrium effects in capillarity and interfacial area in two-phase flow: dynamic pore-network modelling. J. Fluid Mech. 655, 38-71 (2010). For instance, FIG. 1 shows that the capillary pressure-saturation trends converge for different sizes of pore networks with N.sub.B number of pore bodies as the network size increases.

(26) Thus, using REV concepts as shown in Özdemir, M., Özgüç, A., Porosity variation and determination of REV in porous medium of screen meshes. Int. Commun. Heat Mass Transfer 24 (7), 955-964 (1997), it can be stated that Ø.sub.REV=Ø, where Ø.sub.REV is the porosity of REV and Ø is the porosity of the represented porous media. Porosity as used throughout this disclosure is defined as a ratio of void volume to the total volume. Multiplying both sides with respect to the saturation of each container fluid phase leads to Equation 1 below, where V*.sub.i, V* and s.sub.i are the contained fluid volumes within the REV, the total volume of REV, and fluid saturations for each fluid phase hosted within the porous media, respectively.

(27) $\begin{array}{cc}\frac{V_i^{*}}{V^{*}}=\varnothing s_i\mathrm{for}i=n,w.& \mathrm{Equation}1\end{array}$

(28) In various examples, to determine V*.sub.i and V*, REV of a porous media can be represented with a hypothesized homogenous porous medium as shown in FIG. 2. The hypothesized homogenous REV is composed of spherical grains of size r.sub.g, pore bodies of size r.sub.b, and pore throat sizes of r.sub.t as depicted in FIG. 3A. Including such parameters in Equation 1, yields Equation 2 below, wherein N.sub.B, N.sub.t, N.sub.G, and s*.sub.i are the number of pore bodies, pore throats, and grains in the hypothesized homogenous REV, as well as saturations of each of its contained fluid phases. Further, for a hypothesized homogenous REV, it can be assumed in various examples that N.sub.B=x.sub.bN.sub.G and N.sub.t=x.sub.tN.sub.G, which eliminates N.sub.B, N.sub.t, and N.sub.G from Equation 2 as shown in Equation 3 below.

(29) $\begin{array}{cc}\frac{N_B{r}_b^3+N_t{r}_t^3}{N_G{r}_g^3+N_B{r}_b^3+N_t{r}_t^3}=\frac{\varnothing s_i}{s_i^{*}}.& \mathrm{Equation}2\\ \frac{x_b{r}_b^3+x_t{r}_t^3}{r_g^3+x_b{r}_b^3+x_t{r}_t^3}=\frac{\varnothing s_i}{s_i^{*}}.& \mathrm{Equation}3\end{array}$

(30) In the hypothesized REV of FIGS. 2 and 3A, adjacent spherical grains produce inwardly concave pore walls and pore throats. In various examples, this pore throat geometry may be represented using simplified converging-diverging channels as depicted in FIG. 3B. In such examples, the relationship between r.sub.b and r.sub.t can be defined as Equation 4 below based on trigonometry, where β is the orientation angle between pore throat and pore body. Equation 4 may be substituted into Equation 3 to obtain Equation 5 below.

(31) $\begin{array}{cc}{r}_b=\frac{r_t}{1-\tan \left(\beta \right)}.& \mathrm{Equation}4\\ r_t=\sqrt[3]{\frac{\varnothing {s_i\left(1-\tan \beta \right)}^3{r}_g^3}{s_i^{*}\left[x_b+{x_t\left(1-\tan \beta \right)}^3\right]-\varnothing s_i\left[x_b+{x_t\left(1-\tan \beta \right)}^3\right]}}.& \mathrm{Equation}5\end{array}$

(32) Equation 5 relates the properties of the hypothesized REV with the properties of the real porous media. According to the Delaunay Triangulation theorem it is known that N.sub.G≈2N.sub.B, and N.sub.t≈3N.sub.B. Therefore, using Delaunay Tessellation, it can be shown that for a pore network

(33) $x_b=\frac{1}{2}\mathrm{and}{x}_t=\frac{3}{2},$
Equation 5 may be modified to Equation 6 below.

(34) $\begin{array}{cc}{r}_t=\sqrt[3]{\frac{\varnothing {s_i\left(1-\tan \beta \right)}^3{r}_g^3}{s_i^{*}\left[\frac{1}{2}+\frac{3}{2}{\left(1-\tan \beta \right)}^3\right]-\varnothing s_i\left[\frac{1}{2}+\frac{3}{2}{\left(1-\tan \beta \right)}^3\right]}}.& \mathrm{Equation}6\end{array}$

(35) Changes in pore throat size R.sub.t occupied by the fluid-fluid interface as it moves through the complete porous medium, rather than REV, may then be evaluated. In some instances, fluid may flow in drainage conditions in which the non-wetting fluid invades the porous medium and displaces the wetting fluid. During drainage, pore throat size R.sub.t progressively decreases with s.sub.w as shown in FIG. 4. In various examples, the finite difference method can be used to compute R.sub.t in Equation 7 below as the saturation of fluids changes from a range of s.sub.w=s.sub.wi to s.sub.w=1 where s.sub.wi is the irreducible wetting phase saturation and s.sub.n-1 is the saturation of non-wetting fluid at the previous step. The corresponding saturation of fluids can be computed according to Equation 8 and Equation 9 below.

(36) $\begin{array}{cc}{\stackrel{\_}{R}}_t\left(s_n\right)={\stackrel{\_}{R}}_t\left(s_{n-1}\right)+\frac{\delta {\stackrel{\_}{R}}_t}{\delta s_{n-1}}\left[s_n-s_{n-1}\right].& \mathrm{Equation}7\\ s_w=s_{w-1}+\frac{1-s_{\mathrm{wr}}}{t}.& \mathrm{Equation}8\\ s_n=1-s_w.& \mathrm{Equation}9\end{array}$

(37) In Equation 9, t is the number of capillary pressure-saturation data points to be used in generating the capillary pressure-saturation relationship. The estimation of R.sub.t(s.sub.n) depends upon t, and larger values oft may result in more accurate estimates of R.sub.t (s.sub.n). In some examples, to estimate

(38) $\frac{\delta {\stackrel{\_}{R}}_t}{\delta s_{n-1}}$
in Equation 7, it may be assumed that

(39) $\frac{\delta {\stackrel{\_}{R}}_t}{\delta s_n}$
is mathematically similar to

(40) $\frac{\delta r_t}{{\delta }_{s_n}^{*}}$
and accordingly

(41) 0$\frac{\delta {\stackrel{\_}{R}}_t}{\delta s_n}$
may be directly computed from Equation 6, as shown in Equation 10 below. Further, Equation 10 may be incorporated into Equation 7 to obtain Equation 11 below. With regard to Equation 11, in some examples, it may be assumed that s.sub.n*≈s.sub.n for purposes of evaluating the complete porous medium, not solely REV. At s.sub.w=s.sub.wr, R.sub.t(s.sub.n-1) equates to the smallest pore throat within the studied porous medium r.sub.tmin.

(42) $\begin{array}{cc}\frac{\delta {\stackrel{\_}{R}}_t}{\delta s_n}=-\frac{\sqrt[3]{\varnothing {s_n\left(1-\tan \beta \right)}^3{r}_g^3}}{3}\times \left[\frac{\frac{1}{2}+\frac{3}{2}{\left(1-\tan \beta \right)}^3}{{\left[\begin{array}{c}{s}_n^{*}\left[\frac{1}{2}+\frac{3}{2}{\left(1-\tan \beta \right)}^3\right]-\\ \varnothing s_n\left[\frac{1}{2}+\frac{3}{2}{\left(1-\tan \beta \right)}^3\right]\end{array}\right]}^{\frac{4}{3}}}\right].& \mathrm{Equation}10\\ {\stackrel{\_}{R}}_t\left(s_n\right)={\stackrel{\_}{R}}_t\left(s_{n-1}\right)-\frac{\sqrt[3]{\varnothing {s_{n-1}\left(1-\tan \beta \right)}^3{r}_g^3}}{3}\times \hspace{1em}\left[\frac{\frac{1}{2}+\frac{3}{2}{\left(1-\tan \beta \right)}^3}{{\left[\begin{array}{c}{s}_{n-1}\left[\frac{1}{2}+\frac{3}{2}{\left(1-\tan \beta \right)}^3\right]-\\ \varnothing s_{n-1}\left[\frac{1}{2}+\frac{3}{2}{\left(1-\tan \beta \right)}^3\right]\end{array}\right]}^{\frac{4}{3}}}\right]\times \left[s_n-s_{n-1}\right].& \mathrm{Equation}11\end{array}$

(43) In some examples of the present disclosure, variations in pore shape associated with porous media may be ignored. For instance, in such examples, pore bodies and pore throats may be treated as spherical with circular cross-sections. In other examples, a roughness coefficient may be introduced to take into account the highly angular pores, grooves, edges, and tortuous pathways of the porous media that can have a strong impact on the capillary pressure-saturation relationship. The roughness coefficient f.sub.r(s.sub.w) adjusts the pore throat size R.sub.t to obtain an effective pore throat size {circumflex over (R)}.sub.t, which are related according to Equation 12 below. Inclusion of the roughness coefficient may increase the accuracy of the capillary pressure-saturation relationship generation in some instances.
R.sub.t=f.sub.r(s.sub.w)R.sub.t  Equation 12.

(44) From the pore-scale perspective, under drainage conditions, the fluid-fluid interface penetrates into the porous media via piston-like displacement. Joekar Niasar, V., et al., Simulating drainage and imbibition experiments in a high-porosity micromodel using an unstructured pore network model. Water Resour. Res. 45 (2) (2009). Using heuristic evaluation, for piston-like displacements it can be hypothesized that f.sub.r(s.sub.w) varies with s.sub.w in accordance with Equation 13 below, wherein a is a rugosity parameter that will vary depending on the type of porous media being analyzed. For instance, the rugosity parameter may vary depending on the surface roughness of the porous media being analyzed.

(45) $\begin{array}{cc}{f}_r\left(s_w\right)~\frac{1}{\left(s_w^a+{\left(1-s_w\right)}^2\right)\exp \left(s_w\right)}.& \mathrm{Equation}13\end{array}$

(46) FIG. 5 illustrates an example schematic of f.sub.r as a function of s.sub.w under drainage conditions. Initially, as the interface enters the angular capillary its shape is quite regular; however, as the displacement proceeds it will start to distort according to the geometry of the pore, and therefore the roughness coefficient f.sub.r(s.sub.w) increases with a decrease in the saturation of the wetting fluid s.sub.w. As shown in FIG. 5, at maximum f.sub.r(s.sub.w) the interface is similar to the geometry of a capillary. Accordingly, the maximum f.sub.r(s.sub.w) may correspond to the entry capillary pressure. In some instances, f.sub.r(s.sub.w) may be at a maximum when the interface has attained the entry capillary pressure. After reaching the entry capillary pressure, the interface may regain its regular form and break into individual arc meniscus as shown in FIG. 5, meaning that f.sub.r(s.sub.w) declines. FIG. 6 illustrates example trends of f.sub.r(s.sub.w) in response to variations in the rugosity parameter a.

(47) Including Equation 13 in Equation 12 results in Equation 14 below. Including Equation 14 in the Young-Laplace equation generates the capillary pressure-saturation relationship for porous media under drainage conditions, as given by Equation 15 below, where a is the interfacial tension along the fluid-fluid interface, and θ is the contact angle formed between the interface and the pore wall. Data may be input into Equation 15 in order to generate a capillary pressure-saturation relationship for a given porous media under drainage conditions. For instance, Equation 15 may be implemented into a software application or may be utilized with a data analysis tool to generate the pressure-saturation relationship for a given porous media.

(48) $\begin{array}{cc}{\hat{R}}_t=\frac{{\stackrel{\_}{R}}_t\left(s_n\right)}{\left(S_w^a+{\left(1-s_w\right)}^2\right)\exp \left(s_w\right)}.& \mathrm{Equation}14\\ p_c\left(s_w\right)=\frac{2\sigma \cos \left(\theta \right)}{{\hat{R}}_t}.& \mathrm{Equation}15\end{array}$

(49) In some instances of the present disclosure, fluid may flow in imbibition conditions, in which wetting fluid displaces the non-wetting fluid. Generating a capillary pressure-saturation relationship for imbibition conditions is similar to the above-described drainage conditions, with some modifications discussed in the following description. The saturation of fluids may be modified during imbibition conditions as Equations 16 and 17 below, where s.sub.wr is the saturation of wetting fluid at the residual saturation of non-wetting fluid. In various examples, t in the case of imbibition may be similar to that used for drainage displacement.

(50) $\begin{array}{cc}{s}_w=s_{w-1}+\frac{s_{\mathrm{wr}}-s_{\mathrm{wi}}}{t}.& \mathrm{Equation}16\\ s_n=1-s_w.& \mathrm{Equation}17\end{array}$

(51) During imbibition, corner flow is the dominant pore-scale invasion protocol. Therefore, the roughness coefficient f.sub.r(s.sub.w) for the imbibition case may be modified in accordance with Equation 18 below. FIG. 7 illustrates a schematic showing variation in f.sub.r as a result of s.sub.w for the corner flow mechanism. Initially, during corner flow, the individual arc meniscus moves away from the vertex of the angular capillary until they coalesce. Therefore, Jr declines as s.sub.w increases in FIG. 7. After arc interfaces coalesce the probability of the non-wetting phase to get trapped increases. Moreover, it is known from the literature that trapping of non-wetting phase during imbibition is closely associated with the angularity of pores. Accordingly, f.sub.r(s.sub.w) increases with s.sub.w after the individual arc meniscus coalesce. Ma, S., et al., Effect of contact angle on drainage and imbibition in regular polygonal tubes. Colloids Surf. A 117 (3), 273-291(1996); Mason, G., Morrow, N., Capillary behavior of a perfectly wetting liquid in irregular triangular tubes. J. Colloid Interface Sci. 141 (1), 262-274 (1991).

(52) $\begin{array}{cc}{f}_r\left(s_w\right)~\frac{1}{\left(s_w^{\frac{1}{a}}+{\left(1-s_w\right)}^{\frac{1}{2}}\right)\exp \left(1-s_w\right)}.& \mathrm{Equation}18\end{array}$

(53) Including Equation 18 in Equation 12 modifies the effective pore throat size {circumflex over (R)}.sub.t in the imbibition case to Equation 19 below. During imbibition, {circumflex over (R)}.sub.t and a values may remain similar to drainage conditions. Because it is the pore body that controls the displacement dynamics during imbibition rather than pore throat, {circumflex over (R)}.sub.t is converted to a corresponding effective pore body size {circumflex over (R)}.sub.b in accordance with Equation 20 below. Including Equation 20 in the Young-Laplace equation results in the capillary pressure-saturation relationship for imbibition conditions, as given by Equation 21 below. Data may be input into Equation 21 in order to generate a capillary pressure-saturation relationship for a given porous media under imbibition conditions.

(54) $\begin{array}{cc}{\hat{R}}_t=\frac{{\stackrel{\_}{R}}_t\left(s_n\right)}{\left(S_w^{\frac{1}{a}}+{\left(1-s_w\right)}^{\frac{1}{2}}\right)\exp \left(1-s_w\right)}.& \mathrm{Equation}19\\ {\hat{R}}_b=\frac{{\hat{R}}_t}{1-\tan \left(\beta \right)}.& \mathrm{Equation}20\\ p_c\left(s_w\right)=\frac{2\sigma \cos \left(\theta \right)}{{\hat{R}}_b}.& \mathrm{Equation}21\end{array}$

(55) Consistent with the above description, the variables required to generate the capillary pressure-saturation relationship in accordance with Equations 15 and/or 21 are Ø, r.sub.g, r.sub.tmin, a, σ, and θ. In some examples, β may be computed via Equation 4 by replacing r.sub.t with a mean pore throat of studied porous media and r.sub.b with a mean pore body of studied porous media. In other examples, β may be determined from pore-scale images or experiments. In some instances, Ø and/or r.sub.tmin may additionally or alternatively be determined from pore-scale images or experiments. In various examples, σ, the interfacial tension along the fluid-fluid interface, and θ, the contact angle formed between the interface and the pore wall, are measured from pore-scale images or experiments. In other examples σ and θ may be determined by making assumptions to similar situations. Such pore-scale images or experiments, in some instances, may include thin section photomicrographs, e-ray micro CT volume images, confocal microscopy, or other suitable techniques. In instances in which digital images are utilized to determine input parameters, digital images with a higher resolution may tend to produce more accurate capillary pressure-saturation relationships. In some instances, a morphological concept may be applied during image segmentation to achieve more robust results, for example, the morphological concept disclosed by Silin, D., Patzek, T., Pore space morphology analysis using maximal inscribed spheres. Physica A 371 (2), 336-360 (2006).

(56) In some examples, in order to compute r.sub.g, a guess value may initially be used for r.sub.g in Equation 11. Then, in conjunction with Equation 14, the Equation 22 below may be minimized to obtain a value for r.sub.g. With regard to the rugosity parameter a, a guess value may initially be used in some instances. For example, in some instances a value of a=2 may be initially used. Using the computed r.sub.g from Equation 22,

(57) $\frac{\delta {\hat{R}}_t}{\delta s_w}$
may then be calculated tor all saturations in some aspects. In such aspects, if for any saturation value

(58) $\frac{\delta {\hat{R}}_t}{\delta s_w}<0,$
then the effective pore throat size {circumflex over (R)}.sub.t and the grain size of the hypothesized REV r.sub.g may be recalculated with Equations 11, 14, and 22 with smaller values of a until for all saturations

(59) $\frac{\delta {\hat{R}}_t}{\delta s_w}>0.$

(60) 0$\begin{array}{cc}{r}_g=\min \left\{\frac{{\left[r_t-{\hat{R}}_t\left(s_n=0.5\right)\right]}^2}{r_t}\right\}.& \mathrm{Equation}22\end{array}$

(61) FIG. 8 illustrates a box diagram of an example method 800 for generating a capillary pressure-saturation relationship of porous media under drainage conditions, according to an aspect of the present disclosure. Although the examples below are described with reference to the flowchart illustrated in FIG. 8, many other methods of performing the acts associated with FIG. 8 may be used. For example, the order of some of the blocks may be changed, certain blocks may be combined with other blocks, one or more of the blocks may be repeated, and some of the blocks described may be optional.

(62) At step 802, the example method 800 begins. At step 804, s.sub.w and s.sub.n values are determined from s.sub.w=s.sub.wi to s.sub.w=1 with a desired t value using Equations 8 and 9. The saturation values at which the capillary pressures are calculated are accordingly determined. At step 806, values for Ø, β, and r.sub.tmin are determined. In some instances, Ø, β, and r.sub.tmin may be determined from pore-scale images or experiments. In some instances, β may be determined via Equation 4 by replacing r.sub.t with a mean pore throat of studied porous media and r.sub.b with a mean pore body of studied porous media. At step 808, a pore throat size R.sub.t (s.sub.n) and an effective pore throat size {circumflex over (R)}.sub.t are computed using Equations 11 and 14 with an initial guess value for r.sub.g and a=2. In other examples, a may initially begin with a different value. At step 810, the values for r.sub.g and a are updated until for all saturations

(63) $\frac{\delta {\hat{R}}_t}{\delta s_w}\le 0.$
At step 812, parameters are input in order to generate a capillary pressure-saturation relationship for a porous media with Equation 15. At step 814, the example method 800 ends.

(64) FIG. 9 illustrates a box diagram of an example method 900 for generating a capillary pressure-saturation relationship of porous media under imbibition conditions, according to an aspect of the present disclosure. Although the examples below are described with reference to the flowchart illustrated in FIG. 9, many other methods of performing the acts associated with FIG. 9 may be used. For example, the order of some of the blocks may be changed, certain blocks may be combined with other blocks, one or more of the blocks may be repeated, and some of the blocks described may be optional.

(65) At step 902, the example method 900 begins. At step 904, s.sub.w and s.sub.n values are determined from s.sub.w=s.sub.wi to s.sub.w=s.sub.wr with a desired t value using Equations 16 and 17. In some examples, the t value may be similar to the t value used in drainage conditions. At step 906, values for Ø, β, and r.sub.tmin are determined. In some instances, Ø, β, and r.sub.tmin may be determined from pore-scale images or experiments. In some instances, β may be determined via Equation 4 by replacing r.sub.t with a mean pore throat of studied porous media and r.sub.b with a mean pore body of studied porous media. In some instances, the Ø, β, and r.sub.tmin values may be the same or similar to their respective values in drainage conditions. At step 908, a pore throat size R.sub.t(s.sub.n) and an effective pore throat size {circumflex over (R)}.sub.t are computed using Equations 11 and 19 with an initial guess value for r.sub.g and a=2. At step 910, the values for r.sub.g and a are updated until for all saturations

(66) $\frac{\delta {\hat{R}}_t}{\delta s_w}\le 0.$
In some instances, the R.sub.t(s.sub.n), r.sub.g, and a values may be similar to drainage conditions. At step 912, the effective pore throat size {circumflex over (R)}.sub.t is converted to an effective pore body size {circumflex over (R)}.sub.b with Equation 20. At step 914, parameters are input to generate a capillary pressure-saturation relationship for a porous media with Equation 21. At step 916, the example method 900 ends.

(67) FIG. 10 illustrates an example system 100 for forecasting multiphase flow properties of porous media, according to one aspect of the present disclosure. The example system 100 includes different components that are representative of computational processes, routines, and/or algorithms. In some embodiments, the computational processes, routines, and/or algorithms may be specified in one or more instructions stored on a computer readable medium that, when executed by a processor of the system 100, cause the system 100 to perform the operations discussed below. For example, all or part of the computational processes, routines, and/or algorithms may be implemented by the CPU 102 and the memory 104. In other examples, the components of the system 100 may be combined, rearranged, removed, or provided on a separate device or server.

(68) The example system 100 may include a multiphase flow forecaster 106 that is configured to generate (e.g., forecast) multiphase flow properties 110 of a porous media upon receiving input parameters 108. For example, the multiphase flow forecaster 106 may generate a capillary pressure-saturation relationship of the porous media as the multiphase flow properties 110. In various instances, the multiphase flow forecaster 106 is configured to implement the above-described Equations consistent with the above-described method upon receiving the input parameters 108 to generate multiphase flow properties 110 of porous media. For example, the input parameters 108 may be an interfacial tension along an interface between a wetting fluid and a non-wetting fluid, a contact angle between the interface and a pore wall of the porous media, and a pore throat size. The pore throat size may be based on subparameters including a saturation of the wetting fluid, a saturation of the non-wetting fluid, a porosity of the porous media, an orientation angle between a representative pore body size and a representative pore throat size, a representative grain size, and a rugosity of the porous media.

(69) The presently disclosed systems and methods were validated by testing the results from the disclosed systems and methods against experimental data obtained from literature. FIG. 11 illustrates the results of a sensitivity analysis performed for capillary pressure-saturation relationships generated by the provided methods. FIGS. 12 and 13 show a comparison between the provided method (APNA) and mercury injection capillary pressure (MICP) data of various Berea sandstone samples with differing grain textural properties under drainage conditions (FIG. 12) and imbibition conditions (FIG. 13). The MICP data and the input parameters required to generate the capillary pressure-saturation relationships with the presently disclosed method are provided in Table 2 below and were obtained from Churcher, P., French, P., Shaw, J., Schramm, L., Rock properties of Berea sandstone, Baker dolomite, and Indiana limestone, SPE International Symposium on Oilfield Chemistry (1991). The r.sub.tmin value was assumed to be 10×10.sup.−9 m. β.sup.o was computed using Equation 4 as described above, and the r.sub.g and a values were determined according to the above-described method. Overall, FIGS. 12 and 13 demonstrate that the results obtained by the presently disclosed method (APNA—dotted line) were in close agreement with the equivalent experimental data (MICP—star symbols), with the difference between the two being less than an order of magnitude.

(70) TABLE-US-00002 TABLE 2 r.sub.t* r.sub.b* r.sub.tmin r.sub.g s.sub.wi* s.sub.wr* (μm) (μm) β° Ø* (nm a (μm) Berea-1 0.18 0.40 3.70 39.5 41.2 0.19 10 2 38.0 Bata-2 0.16 0.41 6.10 48.0 40.2 0.21 10 2 52.3 Berea-3 0.14 0.49 9.20 54.0 38.9 0.23 10 2 66.7 Orange- 0.25 0.46 18.40 96.0 38.3 0.26 10 2 158.0 Berea

(71) The presently disclosed method was further evaluated using Thomeer's Hyperbola Model, a known capillary pressure model. The Thomeer's empirical equation can be written as Equation 23 below. In Equation 23, p.sub.d represents the displacement pressure (p.sub.c when s.sub.w≈1), V.sub.b is the fractional bulk mercury saturation, V.sub.∞ is the factional bulk mercury saturation at infinite p.sub.c, and F.sub.g is the pore geometric factor. In Themeer's empirical model F.sub.g is illustrative of textural properties of porous media, and thus it is a fitting parameter that reflects the breadth of pore size distribution that is further an indication of sorting of particles in porous media. Large F.sub.g values demonstrate that the porous media is relatively poorly sorted. Through curve fitting, F.sub.g was determined for the presently disclosed method (APNA) and MICP capillary pressure results shown in FIG. 12. While curve fitting, V.sub.∞ was assumed to be equal to the porosity of the respective porous media. The values of F.sub.g for different types of Berea sand stone are shown in Table 3 below.

(72) TABLE-US-00003 TABLE 3 Porous Media APNA MICP Berea-1 0.41 0.43 Berea-2 0.38 0.40 Berea-3 0.35 0.44 Orange-Berea 0.54 0.56

(73) The presently disclosed method was further evaluated against capillary pressure data involving displacement of water by air in sand packs OS-20, OS-30, GB-1, and GB-2, which vary in pore size distribution. The results are shown in FIG. 14. The OS-20 and OS-30 sand pack data was obtained from Schroth, M., et al., Characterization of Miller-Similar silica sands for laboratory hydrologic studies. Soil Sci. Soc. Am. J. 60 (5), 1331 (1996). The GB-1 sand pack data was obtained from Hilpert, M., Miller, C., Pore-morphology-based simulation of drainage in to-tally wetting porous media. Adv. Water Resour. 24 (3-4), 243-255 (2001). The GB-2 sand pack data was obtained from Culligan, K., et al., Interfacial area measurements for unsaturated flow through a porous medium. Water Resour. Res. 40 (12) (2004). The input parameters required to generate the capillary pressure-saturation relationships with the presently disclosed method were also obtained from the directly preceding disclosures and are provided in Table 4 below. β.sup.o was computed using Equation 4 as described above, and the r.sub.g and a values were determined according to the above-described method. Overall, FIG. 14 demonstrates that the results obtained by the presently disclosed method (APNA—dotted line) were in close agreement with the equivalent experimental data (sand pack results—star symbols), with the difference between the two being less than an order of magnitude.

(74) TABLE-US-00004 TABLE 4 r.sub.t* r.sub.b* r.sub.tmin r.sub.g s.sub.wi* (μm) (μm) β° Ø* (μm) a (μm) OS-20 0.07 103 206 26.5 0.34 60.3 1.2 247.0 OS-30 0.05 76.4 153 26.5 0.34 33.9 1.5 181.0 GB-1 0.13 224 448 26.5 0.34 137.0 1.2 577.0 GB-2 0.07 32.8 180 38.6 0.36 20.9 1.2 84.6

(75) The presently disclosed method was additionally compared against capillary pressure results obtained from pore network modeling (PNM) from sandstones imaged using x-ray micro computed tomography. The input variables shown in Table 5 below for the presently disclosed method were directly estimated from the micro-CT images of rock samples and were computed using Avizo Fire 9.2. The micro-CT images of Clashach, Dodington, and Bentheimer sandstones were obtained from “Resources”, Multi-phase flow in porous media (2019), https://www.mfpmresearch.com/resources.html. To compute β, the geometric mean of the pore throat distribution and geometric mean of the pore body distribution were employed with Equation 4 as described in more detail above. Given the limited resolution of x-ray micro computed tomography scans, r.sub.tmin was assumed to be 10×10.sup.−9 m. The r.sub.g and a values were determined according to the above-described method. The comparison results between the provided method (APNA), PNM, and MICP are shown in FIG. 15 for the Clashach, Dodington, and Bentheimer sandstones.

(76) TABLE-US-00005 TABLE 5 r.sub.t* r.sub.b* r.sub.tmin* r.sub.g s.sub.wi* (μm) (μm) β° Ø* (μm) a (μm) Clashach 0 14.0 53.5 35.9 0.14 10 2 252 Dodington 0 17.7 57.9 34.4 0.21 10 2 238 Bentheimer 0 19.6 53.5 32.1 0.19 10 2 272

(77) The results shown in FIG. 15 indicate that the presently disclosed method (APNA) is able to capture the initial capillary pressure-saturation relationship (“p.sub.c(s.sub.w)”) more accurately than PNM and MICP. For instance, there is considerable difference between the results obtained by PNM and the equivalent MICP dataset during the initial phase of displacement, whereas the presently disclosed method (APNA) is able to capture the initial trend more accurately. The improved accuracy of the provided method over conventional PNM approaches for capillary pressure estimation may be due to the analytical nature of the provided method, which makes it free from numerical instabilities while also providing flexibility for users to tune required input parameters to specific cases.

(78) FIG. 16 illustrates graphs of

(79) $\frac{\delta p_c}{\delta s_w}$
for each rock sample of FIG. 15. In many applications, such as reservoir simulation and solving the Buckley-Leverett equation (the classical macroscopic scale immiscible displacement model) it is

(80) $\frac{\delta p_c}{\delta s_w}$
rather than p.sub.c(s.sub.w) that plays a critical role. The results in FIG. 16 indicate that the

(81) $\frac{\delta p_c}{\delta s_w}$
profile is smooth and continuous for MICP, but is highly irregular and discontinuous for PNM that can be ascribed to numerical inaccuracy. In contrast to PNM, the provided method (APNA) is much closer and consistent with MICP data, thus indicating the provided method's greater accuracy as compared to PNM.

(82) Without further elaboration, it is believed that one skilled in the art can use the preceding description to utilize the claimed inventions to their fullest extent. The examples and embodiments disclosed herein are to be construed as merely illustrative and not a limitation of the scope of the present disclosure in any way. It will be apparent to those having skill in the art that changes may be made to the details of the above-described embodiments without departing from the underlying principles discussed. In other words, various modifications and improvements of the embodiments specifically disclosed in the description above are within the scope of the appended claims. For example, any suitable combination of features of the various embodiments described is contemplated.