Stress-Dependent Pore Deformation Effects on Multiphase Flow Properties of Porous Media

Relative permeability and capillary pressure are the governing parameters that characterize multiphase fluid flow in porous media for diverse natural and industrial applications, including surface water infiltration into the ground, CO2 sequestration, and hydrocarbon enhanced recovery. Although the drastic effects of deformation of porous media on single-phase fluid flow have been well established, the stress dependency of flow in multiphase systems is not yet fully explored. Here, stress-dependent relative permeability and capillary pressure are studied in a water-wet carbonate specimen both analytically using fractal and poroelasticity theory and experimentally on the micro-scale and macro-scales by means of X-ray computed micro-tomography and isothermal isotropic triaxial core flooding cell, respectively. Our core flooding program using water/N2 phases shows a systematic decrease in the irreducible water saturation and gas relative permeability in response to an increase in effective stress. Intuitively, a leftward shift of the intersection point of water/gas relative permeability curves is interpreted as an increased affinity of the rock to the gas phase. Using a micro-scale proxy model, we identify a leftward shift in pore size distribution and closure of micro-channels to be responsible for the abovementioned observations. These findings prove the crucial impact of effective stress-induced pore deformation on multiphase flow properties of rock, which are missing from the current characterizations of multiphase flow mechanisms in porous media.

end-point relative permeability (k r−max ), and one more parameter (λ), which defines any curvature in the dependency of relative permeability on the pore size distribution 19 . A further enhancement of the model was achieved by treating the pore space as a bundle of capillary tubes using fractal (D f ) and tortuosity fractal (D t ) dimensional parameters [46][47][48][49][50] . The idea of fractal distribution of capillary tubes was then used to solve classic relative permeability and capillary pressure equations 24 .
Despite an abundance of experimental studies of relative permeability and capillary pressure from both pore-scale and core-scale studies, the impact of pore deformation induced by effective stress changes on multiphase flow properties remains poorly explored 33 . For conditions of increasing effective stress, existing experimental evidences 29,32 reveal contradictory shifts in the relative permeability curves. Recognition of this contradiction has resulted in several authors conducting experiments to tackle the challenges in core-scale experiments 33,34 and pore-scale physical models 44 , independently. However, a comprehensive study of multiphase flow properties (e.g. relative permeability and capillary pressure) in both core-scale and pore-scale under a wide range of effective stress conditions has yet to be undertaken.
In this study, we use X-ray computed micro-tomography to quantify the structure and shape of the pores, together with the pore size distribution and two metrics (fractal dimension and degree of anisotropy), of an initially water-wet Indiana limestone specimen at atmospheric pressure. Then, through a series of core-flooding experiments using water/N 2 pair phases, the stress-dependent porosity, absolute permeability, relative permeability, and capillary pressure of the same specimen are measured under triaxial isotropic effective stress (10 MPa to 30 MPa) and isothermal (40 °C) conditions. Using the measured stress-dependent porosity and pore strain of the core, we reconstruct 3D pore-scale proxy models of the sample at 10, 20, and 30 MPa effective stress conditions. This approach provides a simultaneous investigation of the pore and/or throat shape alteration and pore size distribution indices (e.g. fractal dimension) under variable stress conditions. Additionally, an analytical model is developed based on the fractal theory to reproduce and interpolate the stress-dependent relative permeability and capillary pressure curves. We find that increasing the effective stress results in a leftward shift in pore size distribution, closure of micro-channels, and modification to the flow path tortuosity. These changes are responsible for a dramatic upward shift in capillary pressure curve and a decrease in gas relative permeability. Executing experiments under raising effective stress conditions and constant gas flow rates, we show that increasing applied capillary pressure at the pore-scale due to an increase in gas pressure provides additional energy for the gas phase to invade smaller pore channels, leading to a subsequent decrease in the irreducible water saturation. These findings emphasize the complexity of the physical processes responsible for geomechanical controls on multiphase flow properties in porous media. Based on this research, we discuss the significant role of the stress-dependent relative permeability, capillary pressure, and irreducible wetting phase saturation for some interesting applications in engineering and natural process.

Stress-Dependent Flow Properties from Core-Scale Experiments
A set of single-phase (water) and two-phase (water/N 2 ) core flooding experiments were conducted on a water-wet Indiana limestone core in a high-pressure triaxial cell under a range of isotropic effective confining stress conditions from 0 to 30 MPa and a constant temperature at 40 °C. The equipment 34 allowed for accurate measurement of porosity, absolute permeability, relative permeability, and capillary pressure at several different effective confining stress conditions. First, the confined carbonate core inside the cell was fully saturated with deionized water using the procedure given in the Materials and Methods. Then, the required confining stress inside the cell was applied and the sample porosity measured. At each stress condition, we performed single-phase water flooding tests at several fixed rates (q w = 1-5 mL/min) to obtain stress-dependent absolute permeability. Afterward, the Hassler relative permeability method 51,52 , in which water and N 2 are simultaneously injected directly into the top of the core through two separate lines, was used to evaluate the steady-state pressure drop across the sample for each phase. These measurements were made at several water fractional flow values (f w = 1, 0.8, 0.6, 0.4, 0.2, 0). For the capillary pressure measurement, we employed the modified stationary liquid method 53 , keeping the wetting phase (water) immobile in the core while imposing several steady-state gas flow rates (q g = 1-5 and 8 mL/min). Subtracting the core outlet pressure at steady-state condition from the gas pressure at the inlet provided an estimate of applied dynamic capillary pressure 54 .

Pore-Scale Observations and Conceptual Proxy Modelling
Before running the core-scale experiments, the carbonate specimen was scanned using a stationary X-ray µCT at atmospheric pressure condition. The resulting 3D pore-scale model was synthetically compacted to achieve a proxy structure for the compacted core that was compatible with the imposed effective confining stress conditions. Figure 1 shows the work flow of the process implemented here. More detailed explanations are provided in the Materials and Methods.

Analytical Modelling of Stress-Dependent Multiphase Flow Properties
As noted in 46 , natural porous media are fractals, a fact which has led many researchers to adapt fractal scaling laws to describe the transport and other pore-related properties. Here, the fractal nature of pore sizes was used to develop analytical models for stress-dependent capillary pressure and relative permeability curves, allowing for a better understanding of flow mechanisms on the scale of pores. Based on the fractal scaling laws, the cumulative size-distribution (N(r)) of pores with diameter (r) greater than or equal to r was assumed to follow the fractal scaling law given bellow 55 : where D f and r max are the fractal dimension and maximum pore diameter in porous media, respectively. Based on this definition, the equation for stress-dependent capillary pressure P c (σ′) can be shown to be The terms P e , ε p , λ = 3 − D T − D f and D T in eq. (2) refer to entry capillary pressure, pore strain, stress-dependent pore size distribution index and tortuosity fractal dimension, respectively, and the index i indicates the initial stress condition (σ′ = 0) in this study. Based on the non-linear stress-strain equation for the poroelastic materials, ε p is related to the changes in effective confining stress (σ′) using the following equation 56 , assuming pore strain (ε p ) to be equal to the volumetric strain (ε v ).
where A, B, and C are fitting constants. Haghi et al. 34 defined the fitting constants in eq. (3) as A = (1 − γ s )/K H , B = γ s , and C = K S , where K H and K S are the bulk moduli of hard and soft parts of the rock, respectively. γ s refers to a ratio of the soft part volume to the bulk volume at an unstressed condition. Mixing eq. (2) with Burdine's empirical equation 19 , the following semi-analytical equations are developed for relative permeability curves k r (σ′). where, S wir and k rg−max are irreducible water saturation and maximum gas relative permeability, respectively. A detailed derivation of the above equations is given in SI Appendix A.

Results and Discussions
Core-Scale Approaches. The main contributions of the current research at the core-scale are the experimental results illustrated in Fig. 2 of the stress-dependent structural, single-phase flow, and multiphase flow properties of a porous carbonate rock under a wide range of isotropic effective confining stress conditions. In Fig. 2, circles and squares represent experimental data and colors indicate various effective confining stress conditions. Figure 2(a) shows the decline in stress-dependent permeability (k) versus stress-dependent porosity (ϕ) in response to an increasing effective stress condition, where the colored curve is fitted on the experimental data using the following equation based on the Carman and Kozeny correlation 57 .
Here, τ defines the flow path hydraulic tortuosity inside the porous medium. The black circles in Fig. 2(b) represent the calculated normalized tortuosity (τ τ τ = / i ) from eq. (6) (at known stress-dependent ϕ and k points) and the experimental pore strain data of the sample. The fitted curve is plotted based on the following equation 34 , where the grain shape factor (a) and material constant (m) are calculated using the least squares regression technique ( Fig. 2(b)). This plot clearly reveals an increasing intensity of flow path tortuosity under increasing effective stress conditions. This confirms the simultaneous dual-effect of stress-induced deformation on both pore and connecting-channels sizes, where these are comprehensible through pore strain and hydraulic tortuosity variation, respectively. Fig. 2(c,d) illustrate the core-scale experimental results of how the stress-dependent multiphase flow properties vary under a range of effective confining stress condition from 10 MPa to 30 MPa. These two plots are an experimental demonstration of the remarkable controls that geomechanical processes impose on the multiphase flow properties of a porous carbonate rock. Herein, eq. (2) and (4)(5) are applied as the curve fitting and interpolating equations. More details are given in SI Appendix A.
The results illustrated in Fig. 2(c,d) show leftward (decreasing S wir ) and upward shifts in the relative permeability and capillary pressure curves, respectively. These shifts are conceptually an indication of a decrease in pore sizes, which will be investigated in the next section.
Intuitively, one would expect that a porous medium's increasing affinity to the gas phase should be responsible for the leftward shift of the relative permeability curves and decrease in S wir . However, in our experiments, the www.nature.com/scientificreports www.nature.com/scientificreports/ interfacial tension and fluid properties were constant. To provide quantitative insight into the increasing affinity of the porous media to the gas phase, we provide important evidence from the experimental data at each stress condition: (1) irreducible water saturation ( Fig. 3(a)), (2) water saturation at the intersection point of water and gas relative permeability curves (S wm , Fig. 3(a)), and (3) maximum gas relative permeability ( Fig. 3(b)). These observations, together with our knowledge of the non-changing surface and interfacial tension between the phases and solid during the experiments, reveal the remarkable influence of pore and/or channel deformation on our understanding of the wettability of porous media. Figure 3(c) illustrates an increase in dynamic capillary pressure end-points as effective stress is increased at constant flow rate conditions. The classical capillary number, defined as the ratio of the viscous forces to the capillary forces (i.e. 58 , increases slightly as the core cross-sectional area (A) decreases in response to an increase in effective stress. For capillary number calculations, gas viscosity (μ = . .  58 ; this is consistent with our experimental observations, in which S wir declines as effective stress increases ( Fig. 3(a)). However, given the limited range of calculated Ca for our experiments, in the next section further investigations are warranted for a more representative definition of the capillary number at the pore-scale with more reasonable water capillary desaturation values.
Pore-Scale Approaches. In this section, pore-scale manifestations are evaluated to postulate physical mechanisms supporting the significant core-scale stress-dependent multiphase flow properties measured in our experiments. The proxy modeling procedure given in Fig. 1 is used to develop conceptual 3D models of deformation at the pore-scale.
A qualitative demonstration of stress-induced capillary desaturation of the water phase in a single pore of our carbonate sample is presented in Fig. 4. Here, different colors indicate different effective stress conditions; the color bar corresponding to each is shown on the right side. 3D images of the selected pore at 0 MPa (gray) and 30 MPa (red) effective stress conditions are compared in the same frame in Fig. 4(a). The section in the solid black rectangle is magnified in Fig. 4(b-e) to depict the contraction of a channel connecting the main body to a branch pore with effective stress evolution from 0 MPa to 30 MPa.
Increasing effective stress, a simultaneous increase in both maximum applied capillary pressure (Fig. 3(c)) and channel capillary pressure occurs as the capillary radius r decreases ( Fig. 4(b-e)) following the Young-Laplace equation ) becomes equal to the channel capillary pressure at a critical Ca, the gas viscous pressure (ΔP g ) becomes sufficient to expel, counter-currently, the water phase inside the pore through the channel. This leads to capillary desaturation of the water phase inside the pore. Gas viscous pressure is at its highest level at maximum applied capillary pressure (P cmax Applied ) under each effective stress condition where water saturation decreases to its minimum value at S wir . Scaling down the classical Ca definition for a pore-channel in a porous media, the terms viscous forces and capillary forces can be replaced with viscous pressure and capillary pressure of the channel, respectively 59 . At critical state, maximum gas viscous pressure (ΔP gmax ) at maximum flow rate and steady-state condition, where water phase is immobile, is assumed to be equal to P cmax Applied34,54 .  Fig. 4(f), which provides a multifunctional phase redistribution diagram under variable effective stress conditions, illustrates new insights for stress-dependent critical capillary numbers associated with the pore-channel radius of the carbonate specimen used in this study.
At each effective stress condition, the threshold radius of the medium's channels, in which the capillary desaturation of water starts, can be obtained at = Ca 1 c . Additionally, the transition of a single channel, through its www.nature.com/scientificreports www.nature.com/scientificreports/ deformation path, from a gas capillary barrier to a gas flow conduit (or reverse) at different stress conditions can be recognized. With respect to the channel illustrated in Fig. 4(b-e), increasing effective stress from zero to 10 MPa (Fig. 4(b,c)) results in a transition in the state of the channel from a gas barrier to gas flow conduit; this leads to capillary desaturation of the water phase inside the branch pore. This microscopic transition explains the macroscopic decrease in S wm and S wir in response to an increase in effective stress (Fig. 3(a)). A decrease in S wm and S wir can be translated into an increasing affinity of porous media for the gas phase by increasing effective stress.
To provide further quantitative description of stress-dependent changes in pore-scale structure, pore size distribution, with pore radius (r) being ≥4 μm (Fig. 5(a)), and two metrics, namely 1) 3D fractal dimension (D f ) which quantifies how 3D objects fills pore spaces (Fig. 5(b)) and 2) degree of anisotropy (D A ), which is a measure  www.nature.com/scientificreports www.nature.com/scientificreports/ of 3D symmetry of pores in the media (Fig. 5(c)), are calculated at each effective stress condition. Details on these calculations are given in the Materials and Methods.
As a result of increasing the effective confining stress from 0 MPa to 30 MPa, a leftward shift in effective pore size distribution (PSD) is observed (Fig. 5(a)), while the mean pore radius decreases from 70 μm to 67 μm, respectively. Additionally, the cumulative frequency of pores with r ≤ 52 μm increases from 49% at 0 MPa effective stress to 52% at 30 MPa effective stress; the remaining pores span up to 326.85 μm and 275.25 μm, respectively. These stress-dependent declines in the mean pore radius and size of the largest effective pore, pore connected to the surface of the core, result in an increasing trend of capillary pressure curves and entry capillary pressure, respectively 58 , in response to an increase in effective stress. A similar trend can be realized at the core-scale using the stress-dependent capillary pressure experimental results shown in Fig. 2(d). Additionally, the leftward shift in PSD justifies the decline in maximum gas relative permeability under increasing effective stress conditions ( Fig. 2(b)).
A declining trend was observed in both dimensionless numbers D f and D A in response to an increase in effective stress from 0 MPa to 30 MPa (Fig. 5(b,c)). These two pore structural metrics quantify clearly the rock deformation and its impacts on single-phase and multiphase parameters. The decreasing trend of D f is an indication of more uniform (i.e. well-sorted) pore sizes 58 ; its striking impact on decreasing porosity, absolute permeability, and S wir and increasing slope of capillary pressure curves has been shown experimentally and analytically in the literature 24,[46][47][48][49][50] . In a similar way, our calculated D A represents an increasing state of isotropy of the sample due to pore compaction in response to increasing effective stress.

Implicational Insights
Simultaneous flow of gases and liquids in natural porous media is an indispensable phenomenon to study when we are dealing with subsurface flow processes. Stress-dependent shifts in the relative permeability and capillary pressure curves have a major impact on our classical methods of understanding, quantifying, and modeling multiphase fluid flow in subsurface formations for a wide range of engineering and natural applications, including enhanced oil recovery, CO 2 subsurface storage, heat energy recovery from geothermal reservoirs, and surface water infiltration in the vadose zone.
Although in this study we have conducted the flooding experiments on a standard carbonate specimen using water (wetting phase) and N 2 (non-wetting phase) phases, the discovered systematic stress-dependent behavior of rock could be extended for a wider range of rocks such as sandstone 34 and fluid pairs (e.g. gas + liquid), as gases are typically the non-wetting phase for most natural porous media. For instance, the stress-dependent upward shift in the entry capillary pressure for the studied carbonate rock, at an increasing effective stress condition from 10 MPa to 30 MPa (Fig. 3(c)), almost doubles the injection pressure and energy required to push the gas phase into the porous rock. However, the stress-dependent decrease in irreducible wetting phase saturation resulting from increasing the effective stress condition from 10 MPa to 30 MPa (Fig. 3(a)) leads to an almost 20% increase in the volume of stored gas inside the core and recovery of the wetting phase (i.e. water) from the core. These dramatic changes in gas injection pressure, volume of stored gas, and recovery of the wetting phase are technically and economically significant for similar scenarios such as CO 2 subsurface storage and enhanced recovery of hydrocarbon reservoirs using an miscible/immiscible gas (e.g. CO 2 ) flooding technique. These stress-dependent processes allow us to predict additional CO 2 storage in subsurface formations. In the same way, stress-dependent leftward shifts in the relative permeability curves facilitate rainwater infiltration into the vadose zone and hot water flow through two-phase geothermal reservoirs. This fundamental study paves the way for a more realistic study of fluid flow mechanisms in any natural and industrial applications which are dealing with multiphase fluid flow properties of rock under variable effective stress conditions.

Summary and Conclusions
We have shown experimentally the systematic impact of pore deformation as a result of effective stress changes on the single-phase and multiphase fluid flow properties of a carbonate specimen via a series of core-flooding experiments using water/N 2 phases under isothermal and triaxial isotropic effective stress conditions. Our experiments have revealed leftward and upward shifts in the relative permeability and capillary pressure curves, respectively, with increasing effective stress, and have shown a decreasing shift of S wir , S wm , and k rwmax under increasing effective stress conditions. These experiments have provided us core-scale insights into the linkage between these shifts to decreases in porosity and absolute permeability and increases in pore strain and pore flow tortuosity under an increasing effective stress condition. Fractal and poroelasticity theories were used to drive analytical equations for stress-dependent relative permeability and capillary pressure curves. These equations were used for curve fitting and interpolation in this study.
Using X-ray computed micro-tomography and a proxy modelling technique, we were able to quantify the structure and shape of the pores and channels, pore size distribution, and scaling dimensionless numbers (e.g. fractal dimension and degree of anisotropy) at different effective stress conditions. We have shown that increasing effective stress leads to a leftward shift in PSD and closure of micro-channels. As revealed in the 3D phase redistribution diagram (Fig. 4(f)), we have quantified the stress-dependent threshold radius of the medium's channels at the starting point of capillary desaturation of water in the porous media. Additionally, we have revealed the transition of a single channel, through its deformation path, from a gas capillary barrier into a gas flow conduit in response to an increasing stress conditions, which resulted in the capillary desaturation of the water phase and a macroscopic decrease in S wm and S wir . We have further revealed that leftward shift in PSD is responsible for the increasing trend of capillary pressure curves under an increasing effective stress condition. D f and D A metrics proved quantitatively the role of effective stress condition on pore size distribution and pore isotropy, respectively, in porous media.
These findings underscore the significant impact of effective stress on multiphase flow properties of rock, and manifest the physical mechanisms, such as capillary desaturation and capillary gas blockage due to stress-induced mechanical deformation of pores and/or channels, that control those properties.

Materials and Methods
Core-Scale Experiments. The experiments were conducted on a water-wet Indiana limestone specimen (from Kocurek Industries INC., USA) 3.81 cm in diameter and 10.16 cm in length. The uniaxial compressive strength and Young's modulus of the rock were 25.78 MPa and ≈3 GPa, respectively (Fig. S1). N 2 and deionized water were used as the flooding phases and CO 2 as the displacing phase. Here, an air-water contact angle (α) of 18.8 degrees were measured using a drop shape analyzer (DSA), which indicated the strong water-wet nature of the specimen (Fig. S2). Full re-saturation of the confined sample inside the cell was a key challenge that we tackled using a procedure consisting of 3 steps: 1) flushing the core with 50 pore volume of high pressure CO 2 , 2) vacuuming the sample (at ambient water saturation pressure) from an outlet line over a liquid nitrogen cold trap until a plateau pressure was achieved at the inlet, and 3) injecting high pressure water into the sample and maintaining the pressure for a few hours before re-flushing it with 50 pore volume of water. The experimental procedure in this study is analogous to the protocol described in ref. 34 for Berea sandstone, to which readers are referred for additional details.
Imaging and Proxy Modelling. The X-ray micro-tomography was performed using the micro-CT imaging suite in the Pharmaceutical Orthopedic Research Lab (PORL) at the University of Alberta. The X-ray source within this equipment is a sealed tube with a voltage of 100 kV and a spot size of <5 µm. The imager scans 360 degrees rotational field-of-view at a time using a 12 bit X-ray detector that is fiber-optically coupled to a scintillator. We scanned an 18 mm length of the full-diameter core 5 mm away from the end face of the core. All tomographic images were at 8.6 μm per voxel resolution. The post-processing and quantifications were performed using the CT Analyser program within Bruker microCT 3D Suite software. The proxy modeling process in this study consisted of 7 steps, and the pore stress-strain plot of the rock from the core-scale experiments was an initial input in this cycle. In the first step, we chose an effective stress condition and identified the corresponding pore strain from our core-scale experiments (ε p EXP ). Then the original X-ray µCT images were filtered using median and contrast enhancement techniques. Filtered images were segmented into black and white binary areas representing pores and grains, respectively. Through the thresholding step, a grayscale index smaller than its original value for the rock at zero effective stress was selected manually. Otsu's thresholding method 60 in 3D space was used initially to define the original grayscale index and match the model's 3D porosity with experimentally measured porosity at zero effective stress. The model's pore system was then assigned as a 3D object through bitwise operation, and its volume and strain (ε = − ) were calculated. V P Mi and V P M defined as proxy model's initial and stress-dependent pore volume, respectively. This process was repeated iteratively until the pore strain of the proxy model became equal (±0.1%) to its corresponding experimental value and the representative model used for further 3D analysis.
Pore Size Distribution and Pore Metrics. We implemented the 3D calculation of pore size distribution, fractal dimension, and degree of anisotropy using the CT Analyser program within Bruker microCT 3D Suite Software. For the pore network process, firstly, a "skeletonisation" was performed to identify the medial axis of all pores. Secondly, a "sphere-fitting" measurement was made for all of the voxels along this axis. Pore size is defined as the diameter of the largest sphere which meets two conditions: (1) the sphere encloses the point (but the point is not necessarily the center of the sphere), and (2) the sphere is entirely located within the pore. The Kolmogorov or box-counting method was used to calculate D f in 3D; following this process, the volume was divided into an array of equal cubes and the number of cubes containing part of the object was counted. The process was repeated over a range of box sizes (3 to 100 pixels). Plotting the box sizes (L) versus the number of occupied boxes (N(L)), the slope of log-log linear regression defined D f . Note that a linear relationship between the number of occupied boxes and box sizes in a log-log plot with the slope of D f is also a confirmation of fractal nature of the porous sample based on ∝ − N L L ( ) D f , where this was the case in the current study with the slopes given in Fig. 5(b). Finally, mean intercept length (MIL) analysis was used to measure D A . As part of this process, a grid of lines was sent through the 3D core image containing pores (object) and grains (void) over a large number of 3D angles, and the length of the test line through the analyzed volume was divided by the number of times that the line passed through or intercepted pores. For more details, readers are referred to the software manual (https://www.bruker. com/products/microtomography.html).