Incorporation of Sub‐Resolution Porosity Into Two‐Phase Flow Models With a Multiscale Pore Network for Complex Microporous Rocks

Porous materials, such as carbonate rocks, frequently have pore sizes which span many orders of magnitude. This is a challenge for models that rely on an image of the pore space, since much of the pore space may be unresolved. In this work, sub‐resolution porosity in X‐ray images is characterized using differential imaging which quantifies the difference between a dry scan and 30 wt% potassium iodide brine saturated images. Once characterized, we develop a robust workflow to incorporate the sub‐resolution pore space into a network model using Darcy‐type elements called microlinks. Each grain voxel with sub‐resolution porosity is assigned to the two nearest resolved pores using an automatic dilation algorithm. By including these microlinks with empirical models in flow modeling, we simulate single‐phase and multiphase flow. By fine‐tuning the microlink empirical models, we match permeability, formation factor (the ratio of the resistivity of a rock filled with brine to the resistivity of that brine), and drainage capillary pressure to experimental results. We then show that our model can successfully predict steady‐state relative permeability measurements on a water‐wet Estaillades carbonate sample within the uncertainty of the experiments and modeling. Our approach of incorporating sub‐resolution porosity in two‐phase flow modeling using image‐based multiscale pore network techniques can capture complex pore structures and accurately predict flow behavior in porous materials with a wide range of pore size.

There are two main approaches to multiscale network modeling.In the first, pores and throats at two distinct scales are explicitly incorporated.The smaller elements are captured using ultra-high-resolution imaging (such as nano-CT) and merged with pores and throats imaged, for instance, using micro-CT scanning (Bekri et al., 2005;Jiang et al., 2013;Mehmani & Prodanović, 2014;Mehmani et al., 2013;Prodanović et al., 2015).Generating synthetic pore networks, which explicitly represent pores and throats over multiple scales, has some drawbacks, including potential inaccuracies in representing the actual pore structure, difficulties and uncertainties in fusing synthetic and micro-CT extracted networks, and the possibility of time-consuming and computationally expensive network modeling due to the increased number of elements.
An alternative approach is to replace regions of the pore space where a single image is unable to explicitly resolve pores with Darcy-type elements whose properties are some appropriate average over many smaller elements (Bauer et al., 2011(Bauer et al., , 2012;;Bultreys et al., 2015Bultreys et al., , 2016;;Ruspini et al., 2021;Wang et al., 2022).The use of effective Darcy-like elements that consider the averaged behavior of several combined small-scale elements is more computationally efficient.The problem is how to characterize their properties.One method, used in this paper, is differential imaging which can at least capture locally the connected porosity of unresolved elements (Lin et al., 2016).
Multiscale network models, which represent sub-resolution porosity as a continuous medium or Darcy element, have been used in several studies.Bauer et al. (2011Bauer et al. ( , 2012) ) introduced a Darcy-type throat between two pores only when a macro-throat (a throat that is explicitly resolved in the underlying image) already existed between them.While this successfully captured microporosity as a parallel circuit along the macro-pore throats, it was assumed that connectivity through microporosity only occurs where there is connectivity through resolved macrothroats.This assumption is not valid when the microporous regions provide additional connectivity and accessibility that are not captured by the resolved macro-throats (Lin et al., 2021).Bultreys et al. (2015) extended the approach of Bauer et al. (2012) to include Darcy connections between pores regardless of the presence of macro-throats.They assumed that microlinks exist between pores connecting to the same microporous domain.They tested their model by extracting a multiscale network from an image of Estaillades carbonate and successfully predicted experimentally measured drainage relative permeability from Ott et al. (2015).In a subsequent study, they compared the waterflooding results for Estaillades limestone with experimental results for a Middle Eastern carbonate reservoir rock (Bultreys et al., 2016).Later Wang et al. (2022) used the same approach to create a PNM that effectively handles complex, multiscale porosity.This approach shows promise in addressing sub-resolution porosity using microlinks (connections between explicitly resolved macro-pores).However, they did not validate their multiscale results for waterflooding in Estaillades against experimental data from the same sample.Additionally, their methodology does not guarantee the exclusive assignment of each sub-resolution voxel to a specific microlink.Moreover, there is an artificial constraint on the allowed maximum length of a microlink.In our study, we address these issues by developing a generative algorithm for microlink characterization and proposing a methodology to incorporate wettability in empirical models for microlinks.Then we successfully validate our methodology against a waterflooding experiment on Estaillades limestone.Ruspini et al. (2021) have introduced a multiscale digital rock workflow.They used Darcy-type pores for unresolved porosity.They represented smaller-scale pore structures in a continuum manner using relative permeability and capillary pressure saturation functions calculated from high-resolution images.The workflow was validated using two complex reservoir rocks, and the numerical predictions showed good agreement with experimental data.However, the main concern lies in the acquisition of high-resolution imaging data.While the researchers were able to derive information from high-resolution images for specific areas, this local information may not be representative of the whole sample and might not be easily accessible or feasible in all situations.Finally, wettability is assigned to microporous regions by matching to the measured waterflood capillary pressure which requires additional experimental information.
While existing studies effectively address unresolved porosity in complex porous media, a need persists for a more automated and robust algorithm for characterizing sub-resolution elements and incorporating them into a model with explicitly resolved pores and throats.We have developed a workflow that incorporates sub-resolution porosity from differential imaging into pore network modeling, thus establishing a robust multiscale approach.Unlike previous work, our method incorporates a dilation algorithm to automatically identify microlinks between macro-pores, guaranteeing each accessible microporous voxel is assigned to a specific microlink.This eliminates the need for arbitrary constraints on the microlink length and maintains model efficiency without compromising the accurate representation of the complex pore structure.The following items distinguish the algorithm developed in this paper from previous studies: 1. Based on dilation, we have developed an automatic and robust algorithm that allows us to consider subresolution porosity as microlinks that connect the two nearest resolved pores and construct a multiscale pore network.2. Considering sub-resolution porosity as microlinks, similar to resolved throat elements, means the number of linear equations in the multiscale pore network will be the same as in the resolved network (i.e., the same number of pores).This is in contrast to methods that consider this sub-resolution porosity as a type of pore element, resulting in a larger set of equations to solve for flow calculations.3. Based on our capillary model proposed for all types of wettability, it is possible to derive scanning curves with minimal information regarding the wettability of the system.This eliminates the need for experiments to derive scanning capillary pressure data.
We apply our model to Estaillades carbonate using primary drainage capillary pressure data to calibrate the geometric parameters in the model.We then predict waterflood relative permeability and compare it with experimental data obtained on the same sample.
The structure of this paper is as follows: Section 2 briefly describes the imaging experiment performed on Estaillades carbonate, which serves as a validation experiment.We then elaborate on our workflow for incorporating sub-resolution porosity into pore network modeling.In Section 3, we first assess the sensitivity of the multiscale approach and then validate its results against experimental data.Finally, Section 4 concludes the paper with a summary and discussion of future work.et al. (2019) conducted a study to investigate the behavior of pore occupancy, relative permeability, and flow intermittency measurements using X-ray micro-tomography in a complex carbonate Estaillades sample.The experiment specifically focused on a high capillary number, Ca = 7.3 × 10 6 , which was intentionally chosen to examine the effect of flow intermittency, where the configuration of the phases fluctuates.

Gao
Concurrently, a hitherto unpublished experiment was performed on a different Estaillades sample that had a significantly lower capillary number (2.1 × 10 7 ).In this experiment, the displacement was capillary-dominated with a fixed configuration of phases during steady-state flow.It was conducted on a water-wet Estaillades sample and followed a similar procedure to the high capillary number experiment.
Given its absence of dynamic effects, its water-wet condition, and its notable sub-resolution porosity that has been characterized through differential imaging, this second experiment is ideal for validating our multiscale quasistatic pore network methodology.Subsequent sections will provide a brief overview of the experimental materials, procedures, and analysis.For a more details see Gao et al. (2019) and Lin et al. (2021).

Rock and Fluid Properties
The sample used in this study was Estaillades carbonate.It was cylindrical, with a diameter of 6.00 ± 0.01 mm and a length of 46.70 ± 0.01 mm.A companion sample was composed of approximately 97.9% calcite and 2.1% quartz, as analyzed at Weatherford Laboratories in East Grinstead, UK (Gao et al., 2019).
Using Darcy's law and pressure differential measurements at three flow rates, the absolute permeability of the sample was calculated to be 128 ± 4 mD.The total porosity, estimated from the micro-CT image in Section 2.3, was 28.2 ± 1.5%, with unresolved porosity and macro-porosity accounting for 19.0% and 9.2% respectively.The reported helium porosity for another Estaillades sample was 30.3 ± 0.2%.The image-based porosity value for both this sample and the other Estaillades sample is slightly lower than the helium porosity suggesting that some microporosity might not have been captured.
The non-wetting phase was decane, with a density of 730 kg/m 3 (n-Decane, Acros Organics) and a viscosity of 0.838 mPa s (provided by PubChem, open chemistry database).To achieve optimal contrast between the brine, oil, and rock phases, a 30 wt% potassium iodide (KI) brine was used as the wetting (aqueous or brine) phase.The density of brine was measured to be 1,263 ± 2 kg/m 3 at ambient conditions by weighing a 1 mL drop of the liquid.The viscosity of the brine was determined to be 0.82 ± 0.01 mPa s.
The interfacial tension between brine and decane was measured using the pendant drop method (Andreas et al., 2002;Stauffer, 1965), yielding a value of 47 ± 2 mN/m at ambient conditions, as determined by a Ramé-Hart apparatus (590 F4 series).

Experimental Procedure for In Situ X-Ray Tomography
Our experiment used a 1 mm thick, X-ray transparent Hassler-type flow cell made from carbon fiber epoxy specially designed to withstand high pressures while remaining nearly transparent to X-rays.The sample was positioned within this cell, encased within a Viton sleeve, and connected to the fluid flow-lines via metal fittings.Simultaneous injection of brine and decane was performed through separate ports at a maintained total flow rate of 0.02 mL/min corresponding to a Darcy velocity of 0.71 mm/min.We used this total flow rate and the mean viscosity value of oil and brine to calculate the capillary number, obtaining a value of 2.1 × 10 7 .This categorized our experiment within the capillary dominant regime, common to reservoir settings.The ratio of the Darcy velocity of the aqueous phase to the total velocity of both oil and aqueous phases (f = q w q t ) defines the fractional flow.In the course of the experiment, we increased the fractional flow sequentially to reach steady-state at the following values: 0, 0.15, 0.3, 0.5, 0.7, 0.85, and 1.For further details on the experimental design, refer to Gao et al. (2019).
The experimental procedure was akin to other studies on intermittency by Gao et al. (2019).The steps taken are explained in further detail in Appendix A.
We imaged the fluid configurations using a Zeiss XRM-510 X-Ray microscope.The process employed a flat panel detector, with a voxel size of 3.58 μm.It involved an X-ray energy of 75 keV, an exposure duration of 0.5 s, and 1,601 projections.These projections were converted into a three-dimensional image using the proprietary software of the Versa system.To broaden the vertical (flow direction) field of view for analysis, two images were captured with 25% overlap, stitched together, and cropped into cylindrical forms to remove the sleeve visible at the core boundaries.The final image had dimensions of 1,608 × 1,598 × 2,720 voxels (5.76 mm in diameter and 9.74 mm in length), representing a total volume of 252.1 mm 3 .To ensure a uniform orientation, all the reconstructed images were registered to the dry scan image.The resampling of these images was facilitated by the Lanczos algorithm (Burger & Burge, 2022).
Images were segmented into grains, resolved pores, and sub-resolution pore spaces; the use of 30 wt% KI brine intensified the contrast and emphasized sub-resolution porosity (Lin et al., 2016).Commercial software, Avizo, applied interactive thresholding for segmentation.Figure 1 represents histograms of grayscale CT numbers or intensity values, indicative of X-ray attenuation, for the various phases.These histograms allow for the distinct identification of solid grain, resolved pore, and voxels containing sub-resolution porosity.
Dry and saturated scans, illustrated in Figures 2a and 2b, respectively, leveraged brine's high X-ray adsorption for clear imaging.Differential imaging then enabled the discrimination of micro-and macro-porosity.A non-local means filter was applied to the differential images to reduce noise (Buades et al., 2005(Buades et al., , 2008)), and intensitybased thresholding segmented the image into distinct phases (Lin et al., 2016) (see Figure 2d).
In the context of a dry scan, the average intensity Īdry is related to the porosity (ϕ) and the specific intensities of air (I air ) and rock (I rock ): (1) Water Resources Research For a brine-saturated scan, the average intensity Ībrine is a function of the porosity and the intensities of brine and rock: Porosity is determined by the difference in average intensity between the brine-saturated scan and the dry scan, normalized by the difference in specific intensities of brine and air: This approach essentially measures the contrast change in the image due to the replacement of air by brine in the pore space.
For saturation, the formula is based on the average intensity of a fluid mixture scan Īf w .It is computed by accounting for the porosity, water saturation (S w ), and the specific intensities of brine, oil, and rock:  The water saturation is then calculated from the difference in average intensity between the brine-saturated scan and the fluid mixture scan, normalized by the difference in specific intensities of brine and oil: This process essentially quantifies the contrast change due to the partial replacement of brine by oil.
Importantly, we normalize the images to have the same intensity for rock in brine-saturated, dry, and fluid mixture images and also the same intensity for brine in both brine-saturated and fluid mixture images.This ensures consistent interpretation of image intensities.Taking into account the uncertainty inherent in determining the intensity of phases allows us to estimate the resulting uncertainty in the calculated values of porosity and saturation.

Relative Permeability and Uncertainty Quantification
The pressure difference across the rock sample at steady-state conditions was measured using a differential pressure transducer.We recorded the mean pressure differentials at different fractional flows during the last 2 hours, as well as the corresponding standard deviation.We subtracted the pressure drop in the flow lines themselves (see step 10 of the experimental procedure in Appendix A).The relative permeabilities were calculated by: k rw and k ro represent the relative permeabilities of water and oil, respectively.q w denotes the Darcy velocity of brine (flow rate per unit area) in m/s, while q o denotes the Darcy velocity of oil in m/s.The parameters μ w and μ o signify the viscosities of water and oil respectively, L is the length of the porous medium through which fluid flow occurs, and Δp is the pressure drop across the sample.K is the absolute permeability (m 2 ) determined from singlephase flow (water) through the porous medium based on Darcy's law: where q t is the Darcy velocity for single-phase flow.
The pressure drops in Equations 6 and 7 are corrected by excluding the pressure drop across the tubing in the absence of the sample.To accommodate the impact of local fluctuations in the average saturation in the flow direction, an adjustment is made to the relative permeability according to the methodology proposed by Zhang et al. (2023aZhang et al. ( , 2023b)).

Multiscale Pore Network Modeling Workflow
We assume that any voxel with sub-resolution porosity-whether directly connected to the pore space or through other sub-resolution pores-should be included in the multiscale pore network and considered during the microlink identification process.The workflow is illustrated in Figure 2. Traditional image-based pore network modeling relies on segmented images from dry scans, with extracted networks composed of pores and throats.The first row in Figure 2 illustrates this process, where (a) is the dry scan image, (c) is the segmented dry scan image, and (e) shows the topology of the extracted network.Even for rocks with permeability of 100s of milidarcies some pores and throats may not be fully resolved however (Leu et al., 2014;Saxena et al., 2017Saxena et al., , 2019)).
However, we acquire an additional image using high concentration brine as shown in Figure 2b.Using differential imaging and Figures 2a and 2b, we distinguish voxels with sub-resolution porosity from solid voxels and determine the porosity of each voxel (see Figure 2d).
Those voxels with sub-resolution porosity should be incorporated into the extracted network shown in Figure 2e.This is accomplished using the detailed workflow that will be described later, considering a new element called a microlink.This results in a multiscale network as shown in Figure 2f.

Development of an Automatic Algorithm for Microlink Identification
We treat voxels with sub-resolution porosity as microlinks connecting two resolved pores, provided that these pores are interconnected via these microlinks.We identify the two nearest accessible pores for each voxel with sub-resolution porosity, a task accomplished through an automated dilation algorithm.
Every void space voxel in the dry image is assigned to a pore.We then dilate the image, such that the pore labels start to grow into neighboring solid with sub-resolution voxels, and this process continues until each voxel with sub-resolution porosity has two pore labels assigned.Further details can be found in Appendix B. Figure 3 shows the multiscale network model used in the work, representing Estaillades limestone.

How Displacement Is Modeled
The code to simulate multiphase flow assumes capillary-controlled displacement and has two main steps (Valvatne & Blunt, 2004;Raeini et al., 2018).First, we determine the distribution of fluid phases at a specific pressure for the invading phase.Then, we calculate the pressure or potential field in the pore network by solving the mass or current balance equations for each phase.
Gradually increasing the pressure of the invading phase at the inlet leads to updates in the fluid-fluid interface locations.These updates cause changes in the distribution or saturation of the fluid phases in the network.Displacements of the fluid phase within the center of the element can occur due to potential invasion events, such as piston-like displacement, snap-off, or pore body filling.They can also happen as the fluid interface moves within each individual element, to reach a new equilibrium state.These changes ultimately affect the curvature of the interface, which balances the pressure difference between the phases, where P c is the capillary pressure, P 1 and P 2 are the fluid pressures of the two phases (Phase 1 is brine and Phase 2 is oil), σ is the interfacial tension, and κ is the total curvature of the interface.When the pressure difference exceeds the capillary pressure required for a specific event (piston-like displacement, pore body filling, snap-off) in the pore or throat, the event occurs (Blunt, 2001;Blunt et al., 2002;Valvatne & Blunt, 2004).
For a known distribution of phases, volume conservation is imposed at each pore, assuming incompressible fluids: where j runs over all throats connected to pore i.The flow rate, Q (volume per unit time), between each pair of neighboring pores i and j is determined by multiplying the conductance between the two pores by the pressure difference between them.
where P is pressure, g p is fluid conductance, and L is the length between the pore centers.The conductance represents the inverse of the resistance, which is the sum of the resistance of two pores and the connecting throat.
Conductance is controlled by the geometrical parameters of pores and throats (Blunt, 2001;Blunt et al., 2002;Valvatne & Blunt, 2004), which are determined from the image: where t indicates the connecting throat.The pore lengths, L i and L j are the lengths from the pore-throat interface to the pore center and L t is the length of throat measured as the length between the pore-throat interface i and porethroat interface j.
Equation 9 leads to a system of linear equations for pressure.Solving this system enables us to determine the flow rate for each phase and, subsequently, calculate the permeability for that phase.A similar system of equations can also be used to compute the electrical current through the pore space in response to a potential gradient, invoking conservation of charge.The equations are identical in form, but g represents the electrical conductivity, Q is the current and P the electrical potential.

Characterization of Microlinks
Microlinks have a finite volume, determined from differential imaging.Consider a microlink, denoted m, connecting pores i and j.Consider a voxel with index k with sub-resolution porosity corresponding to this microlink.
The porosity of this voxel is given by ϕ (m) k .Therefore, the porosity of the microlink, denoted as ϕ m , is defined as: where N voxels m represents the number of voxels with sub-resolution porosity associated with microlink m.The total volume of the microlink, denoted as V T m , is equal to the total volume of the corresponding voxels: where dx represents the voxel size.The pore volume of the microlink, V m , is determined by multiplying the microlink porosity by the total volume: The length of the microlink, denoted as l m , is defined as the distance between pores i and j, subtracting their respective pore inscribed radii R i and R j : Water Resources Research where d(i, j) represents the Euclidean distance between the centers of pores i and j.The cross-sectional area of the microlink, denoted as A m , is determined using: Here, β is considered as a tuning parameter, and the same β value is used for all microlinks.We will show later that A m will be used for microlink conductance calculations.To calculate the microlink flow conductance, we also need the permeability of the microlink, denoted as K m .The permeability is found empirically: in this paper, we use the Kozeny-Carman equation (Carman, 1937;Kozeny, 1927), where d g is the average grain diameter which we will use as an adjustable parameter.

Incorporating Microlinks in Flow Simulation
The introduction of microlinks does not alter the total number of pores in the network, ensuring that the size of the linear system of equations remains unchanged.This is a significant advantage of representing subresolution porosity as microlinks since the number of unknown pressures remains constant.We use Equations 9-11 as before but replace the throat length L t with the length l m from Equation 15 and the lengths of pores i and j with their radii.Note that if both a microlink and a macro-throat exist between pores i and j, we separately calculate two Q p,ij values: one for the microlink and one for the macro-throat, and then add them together.
The microlink hydraulic conductance to phase p is defined as follows: where K m is microlink permeability determined from Equation 17, k rp (S w ) is the relative permeability of phase p at brine saturation S w , A m denotes the cross-sectional area of the microlink, Equation 16, and μ p is the dynamic viscosity of phase p.For electrical conductance we use the Archie equation (Archie, 1942), where R w is the brine resistivity, n is the saturation exponent, and F is the formation factor that is the ratio of the resistivity of a rock filled with water (R o ) to the resistivity of water alone (R w ).We assume: where the constant a represents the tortuosity factor, and b is the cementation exponent.
For single-phase flow simulation, we set the water saturation (S w ) to one.From this, we determine conductances and construct the system of equations for flow and electrical conductivity through the network.Solving these systems of equations provides us with the pressure and potential fields, respectively.Next, we calculate the permeability and formation factor of the network.The parameters d g and β, Equations 16 and 17, are tuning parameters used to match the permeability and formation factor obtained from experiments.
In two-phase flow, we must find the brine saturation at each microlink based on the prevailing capillary pressure.This pressure is applied to all accessible resolved elements to determine the curvature of oil/water interfaces, as well as the saturation for all elements including microlinks.To do this, we use the Leverett J-function, which relates capillary pressure to water saturation in a porous medium (Leverett, 1941): Water Resources Research where J is the Leverett J-function, P c is the prevailing capillary pressure, σ is the interfacial tension between the fluids, K m is the microlink permeability, and ϕ m is the microlink porosity, Equation 12.By assuming an empirical model for the J-function based on saturation, we can determine saturation from the prevailing capillary pressure.
For drainage, we use the a power-law model (Brooks, 1965), J* is the entry pressure: if J(S w ) is less than J*, S e is set to 1.The parameter λ represents the power-law exponent, which determines the shape of the capillary pressure-saturation relationship.The normalized saturation of the wetting phase, denoted as S e , is: where S wr is the irreducible water saturation, and S or is the residual oil saturation (for drainage, S or is set to zero).
During drainage, we increase the prevailing capillary pressure, which causes the water saturation to decrease.This continues until we reach a maximum capillary pressure.Each microlink water saturation at this point is represented by S Hyst w .Since microlinks have varying permeability and porosity, each microlink will have a different water saturation S Hyst w ) at this capillary pressure (as given by Equation 21).
Despite the simplicity of the power-law model, it often provides a good approximation for the capillary pressuresaturation relationship, allowing us to determine saturation for each prevailing capillary pressure in microlinks explicitly.However, its application is limited to water-wet systems.Therefore, in this context, we only use it to describe primary drainage.
Afterward, we decrease the capillary pressure to simulate water injection.We employ a more elaborate model that can accommodate different wettabilities and can explicitly determine the saturation for the prevailing capillary pressure (Foroughi et al., 2022).The Leverett J-function based on this capillary pressure model is as follows: with fitting parameters A, B, and C. Equation 24 can be inverted to find the saturation as a function of capillary pressure: After calculating the normalized water saturation, we convert it to water saturation using Equation 23.Instead of using S or and S wr in Equation 23, we use S nr and S * wr to compute the microlink saturation during waterflooding at each prevailing capillary pressure.Next, we will discuss how we obtain these two values.First, by applying the Land trapping model (Land, 1968), we can predict the residual oil saturation, denoted as S nr , at the end of waterflooding for each microlink as follows: where the coefficient C L is determined by: Water Resources Research 10.1029/2023WR036393 FOROUGHI ET AL.
and where we define S or as the maximum residual saturation.At the start of waterflooding, we require the same J (S w ) as determined by both Equations 22 and 24.To achieve this, we introduce a pseudo residual water saturation S * wr defined as follows: where using the waterflooding model (Equation 25) S Hyst e is defined as: Now that we have determined S nr and S * wr , we use these values instead of S or and S wr in Equation 23to compute the microlink saturation during waterflooding at each prevailing capillary pressure.
To examine the behavior of the Leverett J-function for both drainage and waterflooding, refer to the results illustrated in Figure 4a.The Leverett J-function for primary drainage is determined using Equation 22, while for waterflooding, it is found using Equation 24.
After calculating the microlink saturations, we can predict the microlink flow and electrical conductances using Equations 18 and 19, respectively.For flow conductance, empirical models for the relative permeabilities of microlinks are required.As mentioned earlier, we adopt power-law models to describe the relative permeabilities of microlinks.These expressions are often referred to as modified Brooks-Corey relations (Brooks & Corey, 1964), where k max rw and k max ro represent the endpoint relative permeabilities for water and oil, respectively.α w and α o are power law exponents.To account for hysteresis, we employ the Killough model (Killough, 1976).For illustrative purposes, the relative water and oil permeabilities for various S Hyst w values are depicted in Figures 4b and 4c 22, and our proposed model, Equation 24, for waterflooding.We present the results for three saturations at the end of primary drainage based on our proposed model, Equation 24, considering potential variations in saturation following primary drainage.(b) The water relative permeability using the Killough hysteresis model based on power-law relative permeability (k rw ), Equation 30, during drainage and waterflooding.(c) The oil relative permeability using the Killough hysteresis model based on power-law relative permeability (k ro ) Equation 31, during waterflooding.The model considers three initial saturations at the start of waterflooding, denoted as S Hyst w .

Bypassing Zero Saturation in Two-Phase Conductivity
We encountered a situation where the saturation of the wetting phase in a pore became very low, effectively halting the flow (see Figure 5).To ensure continuity of flow when water is present in microlinks, we devised a bypassing mechanism.This mechanism works by rerouting the flow around the pore with zero or very low saturation for the wetting phase, allowing the flow to continue despite the blockage.To facilitate this bypassing mechanism, we set L ij /g p,ij to L m /g p,ij whenever L ij /g p,ij , calculated using Equation 11, is less than L m /g p,ij .This adjustment allows us to maintain a non-zero conductance.

Results and Discussion
In this section, we will evaluate the performance of our multiscale PNM.The benchmark for this assessment is the Estaillades rock sample, described in Section 2.3.
Estaillades rock is a heterogeneous material characterized by a broad range of pore throat sizes.Figure 6 depicts the throat size distribution obtained from MICP porosimetry, as reported by Tanino and Blunt (2012).This distribution spans from a few nanometers to tens of microns, indicating a wide range of radii.
We use the voxel size as a threshold (red dashed line in Figure 6) to distinguish between resolved and unresolved pore space noting that more than one voxel is required to accurately resolve the width of a pore or throat.For heterogeneous rocks such as Estaillades, the throat size distribution typically exhibits a bimodal pattern: the peak at smaller radii corresponds to micropores and micro-throats, while the peak at larger radii represents macro-pores and macro-throats.By scrutinizing the voxel size threshold and the distribution, we can differentiate between two types of unresolved pore space: micropores associated with the first peak, and macro-pores linked to the second.
Considering the considerable fraction of unresolved macro-pores (see Figure 6), our multiscale modeling approach needs to accurately differentiate them from micropores.
The subsequent Section 3.1 presents single-phase flow results based on this multiscale pore network.Here, we tune the parameters of the microlinks to match the experimental permeability.Next, in Section 3.2, we focus on two-phase flow simulations.Initially, we match the reported primary drainage capillary pressure for Estaillades.Following this, in Section 3.3, we present the drainage relative permeability results.In the end, we display the predicted results for waterflooding relative permeability, comparing them against the experimentally measured values.
The equations to be used in the following sections are summarized in Table 1.

Single-Phase Flow -Model Tuning
In this section, we present the results of single-phase flow analysis conducted on three sub-volumes of the sample.Each sub-volume consists of 1,127 × 1,127 × 1,127 voxels, representing a substantial portion of the rock sample.The analysis focused on determining the permeability and formation factor, which were obtained through simulations.Permeability was calculated using Equation 7 through the simulation of flow within the extracted pore network.The formation factor was determined by simulating the electrical current through the rock and calculating its electrical resistivity (the inverse of electrical conductivity).
To calculate the conductance of microlink, two distinct grain diameters (d g ) have been used: 14.0 and 1.75 μm.This choice is essential to capture two types of unresolved porosity.The selection of these two values allows us to match the measured permeability and primary drainage capillary pressure.For electrical simulations, we use an empirical conductivity model for microlinks based on Archie's law (see Table 1).
The results for single-phase flow are presented in Table 2. To assess the influence of sub-resolution porosity on flow behavior, two simulation types were compared.The first disregarded the presence of sub-resolution porosity, while the second, based on our methodology, considered the microporous nature of the rock.All three subvolumes exhibited similar fractions of resolved and unresolved pore space.When sub-resolution porosity was not considered, the predicted permeability using pore network modeling was consistently one order of magnitude lower than the experimentally measured value.Similarly, the predicted formation factors exhibited were unreasonably high (see Table 2).
However, by incorporating sub-resolution porosity and implementing our multiscale PNM, we successfully calibrated the microlink parameters to obtain reasonable permeability values that aligned with the experimental measurements (see Table 2 for all sub-volumes).Furthermore, the predicted formation factors experienced a significant improvement and fell within the reported range for carbonate rocks, namely 10-100 (Youssef et al., 2008).
The tuning parameters in this study included the constant d g , which corresponds to the grain diameter in the Kozeny-Carman relation, and the parameter β which scales the microlink cross-sectional area in Equation 16.By carefully adjusting these parameters, we managed to achieve more accurate simulations that matched experimental values for permeability: β = 2.5 for all sub-volumes.The tuned value for d g will be discussed in the next section.

Two-Phase Flow
As mentioned earlier, the throat size distribution for Estaillades displays two distinct peaks, corresponding to two types of porosity.Additionally, based on Figure 6 and considering the measured voxel size of the sample (3.58 μm), it is evident that a significant fraction of the pore space in this rock remains unresolved.We assign the grain diameter to each microlink as follows: Here, ϕ m represents the mean porosity of microlinks, and σ ϕ m is the standard deviation of microlink porosity.Additionally, ϕ m,c denotes the porosity criterion value corresponding to each microlink, which is used to assign d g to that microlink.
The N symbol represents a random selection of ϕ m,c from a normal distribution with the mean and standard deviation indicated.Then the d g and hence permeability, Equation 17, of the microlink is determined based on the  32.This probabilistic model allows the capillary pressure to vary smoothly in the transition from resolved macro-pores, to unresolved but larger pores to true microporosity, to reproduce the behavior observed in Figure 7.It is assumed that the pores and throats extracted from the resolved pore space are strongly water-wet, with a contact angle of 0°for drainage.The Leverett J-function used for primary drainage is presented in Table 1.
To validate the performance of our approach, we generated capillary pressure for three sub-volumes of the Estaillades rock sample (Figure 7).We then compared these results with experimental data from two literature sources: Tanino and Blunt (2012) and Alyafei and Blunt (2016).There is excellent agreement between the simulated and experimental results.The simulated capillary pressure values were slightly higher than the data of Alyafei and Blunt (2016), which may be due to the lower permeability observed in our sample which is not fully accounted for by J-function scaling.Alyafei and Blunt (2016) reported a permeability of 182 mD, while the sample studied by Tanino and Blunt (2012) had a permeability of 166 mD compared to our value of 128 mD.

Prediction and Validation
In Figure 8, we depict the primary drainage relative permeability for the three sub-volumes of the Estaillades sample.The graphs in Figures 8a and 8b show the same data on Cartesian and logarithmic y-axes respectively.The relative permeability trends of water and oil are generally similar across all sub-volumes.Here, however, we do not have experimental data for comparison.However, similar results for Estaillades rocks in drainage have been reported by other modeling studies (Bultreys et al., 2015;Carrillo et al., 2022).
The relative permeability model for microlinks (Table 1) assumes water-wet behavior during waterflooding.Pores and throats extracted from the resolved pore space are considered water-wet, with a contact angle of 45°for waterflooding.This contact angle is similar to the one reported for a water-wet Bentheimer sample in our previous studies (Blunt et al., 2019).
The final results of relative permeability for waterflooding in the Estaillades sample are plotted against experimentally determined values in Figure 9.For the sake of clarity, we have presented each subvolume using separate plots.The various lines on the graph represent the inherent uncertainties associated with pore network modeling and different assumptions about initial saturation.For each subvolume, we generated 310 realizations by varying the initial brine saturation and using different seeds for the waterflooding simulations.The initial brine saturation is consistent and within the range of the experiments.The results are shown on both Cartesian (Figures 9a, 9c, and  9e) and logarithmic scales (Figures 9b, 9d, and 9f).The model predictions match the experimental results within the uncertainty bounds of the experiment and simulations.For comparison Carrillo et al. (2022) performed DNS using a multiphase micro-continuum model on images of Estaillades and predicted relative permeabilities for drainage and imbibition.While their results are similar to ours at high water saturations, they predicted a water relative relative permeability of approximately 0.05-0.1 for water saturations around 0.2.In our model, and the experiments, k rw is much lower for both drainage and imbibition, associated with the low flow conductance in wetting layers and microlinks.

Conclusions and Future Work
We have presented a multiscale PNM to simulate fluid flow in heterogeneous rocks with unresolved porosity.The model incorporates Darcy-type elements, microlinks, enabling accurate representation of the pore space and flow behavior.Incorporation of this unresolved pore space in modeling is crucial for obtaining accurate pore volume measurements and ensuring correct accessibility and flow.
We treat sub-resolution porosity as microlinks that connect the nearest resolved pores found using a dilationbased algorithm.This algorithm assigns each grain voxel, which contains sub-resolution porosity, to its two closest resolved pores.Voxel groups with the same two nearest pores are categorized as a microlink connecting them.We employ empirical models to calculate capillary pressure and relative permeability through the microlinks.
We test our model on Estaillades carbonate which has a wide range of pore size.We use differential imaging to identify both the resolved macro-pore space and unresolved porosity.We first tune the model parameters to provide a good match to the measured permeability and drainage capillary pressure.We assign two types of microlinks based on their porosity: the more porous links are assumed to be principally associated with pores whose size is just below the image resolution, while lower porosity elements represent smaller microporosity.
We predict both drainage and waterflooding relative permeabilities for water-wet conditions.The predicted waterflood relative permeabilities fall within the uncertainty of the experimental measurements.
In the future, our intention is to expand the application of this methodology, coupled with our wettability optimization workflow (Foroughi et al., 2020(Foroughi et al., , 2021)), across a diverse range of data sets.
The dilation process is performed randomly across six directions to prevent directional bias.This process only applies to voxels with sub-resolution porosity.Grain voxels without sub-resolution porosity are not part of the dilation process.This procedure is repeated, updating the microPore1 matrix through subsequent dilation levels until no more changes occur.By this stage, all grain voxels containing sub-resolution porosity that are directly accessible to void space or through other microporous voxels will have been assigned the nearest pore label.The output of this process for cases in which we have two and three pores and a gray area (the region with sub-resolution porosity) in between is shown in the second columns of Figures B2a and B2b, respectively.
• Assigning Next-Nearest Pore Labels to Voxels with Sub-Resolution Porosity-An Algorithmic Approach: The second step involves assigning the next-nearest pore label to each grain with sub-resolution porosity.We use a similar algorithm as before.We initialize a new matrix, microPore2, of the same size as microPore1, with all elements initially set to zero.We then dilate microPore1 and assign the results to microPore2.Dilation occurs through void voxels and grains with sub-resolution porosity.The dilation of microPore2 continues through the void voxels and grains with sub-resolution porosity until no further changes occur in microPore2.
After these two steps, we obtain two matrices, each the same size as the dry scan image, in which each voxel with sub-resolution porosity is associated with its nearest and next-nearest pore, as indicated by microPore1 and microPore2, respectively.A schematic representation of next-nearest pores for cases where we have two or three pores with a gray area (the region with sub-resolution porosity) in between is shown in the third columns of Figures B2a and B2b, respectively.• Identification of Microlinks: Once we have successfully assigned two labels to each grain voxel with subresolution porosity-either in direct contact with a void voxel or connected through other voxels with subresolution porosity-we can identify the microlinks.A microlink comprises all grain voxels with subresolution porosity that share the same pair of labels, regardless of their order.In other words, all elements in the microPore1 matrix with label i, and corresponding elements in the microPore2 matrix with label j (or vice versa), correspond to a microlink m that connects Pore i and Pore j.A schematic representation of a microlink between adjacent pores, determined after superimposing microPore1 and microPore2 on each other for cases where we have two or three pores with a gray area (the region with sub-resolution porosity) in between, is shown in the fourth columns of Figures B2a and B2b, respectively.We then assign the parameters of microlink m based on these grouped voxels.Note that many microporous voxels will be assigned to the same microlinks: the number of microlinks is therefore much smaller than the number of voxels in the image.

B1. Addressing Computational Challenges in the Proposed Algorithm
The dilation algorithm poses substantial computational challenges, particularly when processing large images.To mitigate this issue, we implemented a strategy that subdivides the label matrix into smaller sub-matrices.The computations for each of these sub-matrices are executed concurrently using multi-threading, significantly enhancing the efficiency of the algorithm.By employing this strategy, we successfully navigated the computational complexity associated with the dilation algorithm, enabling the processing of large images without encountering excessive computational costs.

Figure 1 .
Figure 1.Histograms of intensity for different phases obtained from the brine-saturated scan.Each line represents a specific phase found in the scan, including the "All phases," category.

Figure 2 .
Figure 2. (a) Displays a 2D slice from a dry-scan 3D image used to identify resolved pore space.The different colors represent the density of each phase.In the dry image, red represents the highest density (grains), green represents the intermediate phase (microporous regions), while blue represents the resolvable pore space.(b) Presents the same 2D slice as in panel (a), but of the rock saturated with 30% potassium iodide brine.This image is used alongside (a) for differential imaging to discern solid voxels containing sub-resolution porosity.(c) Features the same slice as (a) after watershed segmentation, highlighting resolved pores in black for conventional pore network extraction.(d) Depicts the same slice as in panels (a), (b), and (c), but segmented to label three components: resolved pore voxels (black), solid voxels (blue), and voxels containing sub-resolution porosity (gray), which are distinguished through differential imaging.The porosity of the red voxels is determined using the grayscale values from (a) and (b).(e) Presents the network, composed of macro-pores and macro-throats, extracted from the resolved pores (c).(f)Shows the multiscale pore network, which incorporates the extracted network from (e), the pore label image generated by network extraction (not shown here), and the 3D image from (d).This facilitates the identification and characterization of microlinks, resulting in a network composed of macro-pores, macro-throats, and additional microlinks.

Figure 3 .
Figure 3. (a) Topology of the extracted pore network based on the resolved pore space.(b) Updated topology of the same pore network after incorporation of microlinks (referred to as the multiscale pore network).The colors of the pores (spheres) and throats (cylinders) represent their inscribed radii.The microlinks are shown in gray.

Figure 4 .
Figure 4. (a) The Leverett J-function employed in our model.The diagram depicts the power-law model for drainage, Equation22, and our proposed model, Equation24, for waterflooding.We present the results for three saturations at the end of primary drainage based on our proposed model, Equation24, considering potential variations in saturation following primary drainage.(b) The water relative permeability using the Killough hysteresis model based on power-law relative permeability (k rw ), Equation30, during drainage and waterflooding.(c) The oil relative permeability using the Killough hysteresis model based on power-law relative permeability (k ro ) Equation 31, during waterflooding.The model considers three initial saturations at the start of waterflooding, denoted as S Hyst w .

Figure 5 .
Figure5.Illustration of a pore, denoted as j, located between two microlinks (m ij and m jk ), while all other elements belong to resolved porosity.This pore tends to fill before the microlinks (with blue representing the wetting phase, water, and red representing the non-wetting phase, oil), requiring a bypass of pore j to ensure that the flow is directed exclusively through the microlinks, thereby avoiding fluid trapping.

Figure 6 .
Figure 6.The throat radius distribution for the Estaillades limestone is depicted here, highlighting the clear presence of a bimodal pore size distribution.The dashed line represents the image voxel size of 3.58 μm.Data from Tanino and Blunt (2012).

Figure 7 .
Figure 7. Dimensionless measured capillary pressure (written as a dimensionless J function, Equation21) as a function of water saturation and the results of multiscale modeling for three sub-volumes.Data fromTanino and Blunt (2012) andAlyafei and Blunt (2016).

Figure 8 .
Figure 8. Comparative analysis of relative permeability for water and oil across three Estaillades sub-volumes, for primary drainage.

Figure 9 .
Figure 9. Comparative analysis of relative permeability computed on three Estaillades sub-volumes, for waterflooding using our multiscale pore network model, versus experimental results.A good agreement between the modeling and experimental data is evident.The error bars represent the uncertainty in the experiments, taking into account variations in the measured saturation profiles, and measurement uncertainty (see Zhang et al. (2023b) for further details).The blue dotted line corresponds to the experimental relative permeability of water, and the red dashed line corresponds to the experimental relative permeability of oil.Panels (a) and (b) correspond to sub-volume 1, panels (c) and (d) correspond to sub-volume 2, and panels (e) and (f) correspond to sub-volume 3.These modeling results are derived from 310 realizations.

Figure B2 .
Figure B2.The first column presents a schematic of resolved pores with a gray area denoting the voxels that connect these pores via sub-resolution porosity.The first row consists of two pores, while the second row contains three.The identification of these pores is facilitated by a pore network extraction on the resolved pore space.The second column highlights the use of the dilation algorithm, which assists in recognizing the nearest pore label for each voxel with sub-resolution porosity.With further dilation, we can identify the next-nearest pore label for each voxel with sub-resolution porosity, as illustrated in the third column.The fourth column demonstrates the process of superimposing the nearest and next-nearest labels, enabling the identification of microlinks.Each voxel sharing the same set of nearest and next-nearest labels, irrespective of order, corresponds to a unique microlink.

Table 1
Summary of Parameters and Equations Used for Microlinks Is the porosity, K the permeability and F is the formation factor.