1. Introduction
We typically produce only 40% of original oil in place (OOIP) (Gupta et al. [
1]). Thus, understanding the fluid flow mechanisms in the reservoir and what can be done to improve recovery continues to be a critical area of research. Because the fundamental factors that dominate the recovery remain the same, the reservoir models/simulators used to evaluate such recoveries need to be improved so that we can better understand the interaction mechanism of fluid-rock and fluid flow behavior. Usually, core flooding and centrifugal methods are used to determine the capillary pressures, relative permeabilities, recovery factors and the overall recoveries of the flooding processes. Although important to conduct, these experiments are time consuming, limited to a few flooding cycles and require preserved cores from the target formation. Numerically simulated cores mimic the rock and reservoir conditions and thus allow for numerous and rapid simulation of flooding scenarios. However, the accuracy of simulations depends on how well the models represent the actual reservoir.
Relative permeability (kr) is a key property that governs multiphase flow in porous medium and plays a central role in enhanced hydrocarbon recovery processes. The flow of a phase is directly determined by their relative permeabilities in the subsurface reservoir, which impacts the hydrocarbon production rate and ultimate recovery. Accurate relative permeability estimation is critical for improved reservoir characterization as it can affect the outcomes of breakthrough time, history matching, production forecast, and estimated ultimate recovery.
Experimental and modeling studies have tried to understand the controlling factors of relative permeability. Morgon and Gordon [
2] showed that relative permeability was strongly dependent on fluid saturation and pore structure only. Naar and Henderson [
3] studied the impact of saturation history on relative permeability. Fluid properties that also effect the relative permeability include temperature (Sedaee Sola et al. [
4]; Hamouda et al. [
5]), flow rate (Yeganeh et al. [
6]), viscosity ratio (Peters and Khataniar [
7]) and interfacial tension (Cinar et al. [
8]). Empirical models, such as those of Brooks-Corey [
9] and Land [
10], have also been proposed for accurate estimate of relative permeability. However, all these models either account for the effect of only one parameter or require fitting parameters which have limited predictive capabilities.
Khorsandi et al. [
11] recently proposed an equation-of-state approach (EoS) to provide a single-valued estimate of relative permeability based on six state parameters: phase saturation, phase connectivity, fluid–fluid interfacial area, pore structure, wettability of the medium, and capillary number. Mcclure et al. [
12,
13] showed that including phase saturation, connectivity, and fluid-fluid interfacial area could capture hysteresis in relative permeability well. Purswani et al. [
14] then developed a simplified workflow for the development of an EoS for relative permeability as a function of saturation (
S) and phase connectivity (
X) only, by defining relevant boundary conditions and constraining the EoS under limiting conditions. They assumed other state variables such as contact angle were constant. However, their PNM data show differences in relative permeability values for different directional paths that intersect at the same
S and
X values. This difference indicates that saturation and connectivity may not be sufficient to uniquely model relative permeability. Thus, interfacial area might be one of the factors that can reduce the hysteresis and it is worthwhile to include it in the functional form for relative permeability.
This study investigates the sensitivities of the three state parameters (
S,
X and
A) on relative permeability estimation. We also analytically calculate the fluid-–fluid interfacial area (
A) for secondary drainage and add it to the existing dataset (primary drainage and imbibition) in Mukherjee et al. [
15]. These findings can have direct implication on static and dynamic reservoir simulations and the knowledge of the sensitivities in
S-X-A space can help in the formulation of a more accurate and robust state function for relative permeability.
2. Methods
Relative permeability, phase saturation, phase connectivity, and fluid–fluid interfacial area for primary drainage and imbibition were generated using the classical pore network model and analytical technique (Mukherjee et al. [
15]). We have improved on the previous model to generate data for the secondary drainage for this study and further evaluate the sensitivity of the state parameters (
S,
X and
A) on relative permeability estimation. In this section, we present a summary of the key features of the pore network model and the method used to calculate the interfacial area.
2.1. Pore-Network Model
While the experimental techniques to measure the petrophysical properties help constrain the physics during two phase flow in porous media, they cannot be used to explore a complete functional dependency of relative permeability. Pore network modeling (PNM) simulations on the other hand simplify the real systems to run numerous simulations in a short span of time. The network model used in this study (Valvatne et al. [
16]) is based on the pore filling mechanisms, such as percolation or invasion–percolation. These mechanisms govern the flow of the invasive fluid into the pores and consequently affect the capillary pressures. During flow, when an entry pressure to a pore is satisfied, invasion takes place. By tracking the saturation changes for each pore filling event, the relative permeability for each phase is then calculated as the ratio of flow rate of that phase to the total flow rate for each simulation.
Table 1 lists the properties of the network model created using a micro-CT image of a Bentheimer sandstone. Simulations of the model require additional inputs such as interfacial tension (
), viscosity (
, density (
, and wettability (θ) as shown in
Table 2. A total of 200 PNM simulations for primary drainage (PD), imbibition, and secondary drainage (SD) were run at equally spaced initial water saturations from 0.0 to 0.9 along the PD curve at a constant contact angle of 0° for the wetting phase. Small saturation jumps of 0.002 units were used for each simulation run.
2.2. Fluid-Fluid Interfacial Area
From specific outputs from the PNM simulations such as relative permeability (kr), saturation (S), capillary pressure (Pc), location of fluid phases, corner half angles (β) and dimensions of pores and throat elements, fluid-fluid interfacial area is calculated. We have calculated both the corner areas or Arc Menisci (AMs) and throat areas or Main Terminal Menisci (MTMs) for three different types of pore shapes (cube, triangular, and spherical). The total specific fluid–fluid interfacial area (A) is the summation of both the corner areas and throat areas normalized to the bulk volume of the network model.
2.3. Predictor Importance
The out-of-bag predictor importance method was used to determine the relative importance of
S,
X and
A (predictors) on the relative permeability estimates (response). The influence of a predictor increases with its importance, which is represented as an estimated number. If the predictor estimates are high, then permuting its values should affect the model more than permuting the values of a predictor with a lower estimate. The flowing steps describe the evaluation of predictor importance estimates by permutation (
Figure 1):
Creating a predictive model: The random forest ensemble method is used to construct multiple decision trees during training, using all the available data as input (Grover [
17]). Random features are chosen at each step of the decision tree that helps in reducing the variance and improving the accuracy of the model.
Out-of-bag error: For each of the predictor variables, the observations are permuted randomly, and the model error is estimated for each tree. The final estimates are then evaluated by averaging the outcomes of all the decision trees. The benefit of using the random forest ensemble is that it handles overfitting and reduces the variance.
3. Results
In this section, we first discuss the results from the PNM simulations for the S-X paths for a completely water wet rock followed by the specific area curves for the respective simulations. Lastly, the estimates regarding the sensitivity analysis of saturation S, connectivity X, and specific area A on relative permeability are presented.
Figure 2 shows the
S-X plots for five primary drainage, imbibition, and secondary drainage scans for both wetting and non-wetting phases. Each curve starts with the primary drainage (black solid line) to terminate at a specified termination point, followed by the imbibition cycle (blue dashes) up to a residual saturation point. The secondary drainage cycle (red dots) starts from the residual saturations and goes to 100% saturation for the non-wetting phase (0% saturation for the wetting phase).
Figure 3 shows an example of the total specific interfacial area (
A) as a function of the non-wetting phase saturation for a primary drainage (PD), imbibition and secondary drainage (SD) cycle. As we start increasing the non-wetting phase saturation with PD, we observe that for an increase in saturation of the invading phase, the interfacial area increases from zero at
Snw = 0.0, reaches a maximum peak around ~70%
Snw, and then decreases to zero at
Snw = 1.0. The increase in area results from an increasing volume of the invading non-wetting phase in the model. After the peak, the area must decrease to 0.0 as the nonwetting-phase saturation approaches 1.0 because the wetting phase vanishes there.
The imbibition scan starts at the same
A as for PD termination, but as wetting phase imbibes the interfacial area deviates from the PD curve, showing a substantially larger peak. This deviation from the PD curve and the increase in the peak level is the result of nonwetting phase trapping. The scans terminate at the residual oil saturation where relative permeability is numerically close to zero. The area for SD follows the same trend. However, during SD, the saturation of the non-wetting phase increases significantly so that the trapped non-wetting phase reconnects, decreasing the area until it reaches 0.0 at
Snw = 1.0. Ideally the SD curve should start at the termination point of the imbibition cycle. However, in
Figure 3, the SD curves (red) start at a higher saturation than the saturation at which the imbibition curves (blue) terminate. We believe this small discrepancy of an extended tail of the imbibition curve at lower non-wetting phase saturations is an error due to the not-so-accurate values of capillary pressures from the PNM simulations at near-residual saturations.
Sensitivity Analysis
In this section, we discuss the relative importance of the state parameters
S,
X and
A with respect to
kr estimations. The values of the estimates as an output from the out-of-bag predictor importance method are only comparable when we consider the same input data set. For different input datasets, the values of the estimates can vary and thus are not comparable from case to case, although the trends (increasing or decreasing) can still be compared. We input all the data into the model without splitting it into train/test data and
Figure 4 shows a perfect fit of the data by the random forest ensemble. Like Elsheikh et al. [
19], the ensemble also achieves a very high R
2 value of unity.
Figure 5 shows the estimates for the predictor variables
S,
X and
A for the entire data set pertaining to the primary drainage, imbibition, and secondary drainage cycles for both the wetting and non-wetting phases. We observe from the results that saturation is the most important parameter for the wetting phase followed by connectivity and then specific area. For the non-wetting phase, connectivity is the principal factor followed by saturation and specific area. Physically, since the rock studied here is water wet, the wetting phase flows as layers adjacent to the walls of the rock, and thus the amount or phase saturation controls the flow. However, for the non-wetting phase, connectivity plays a critical role due to the snapping off events during imbibition and the reconnection of the isolated phases during secondary drainage.
Further, we can zoom-in on each of the wetting and non-wetting phase data to analyze the sensitivities of the predictors across the saturation domain. The
S-X space is divided into four quadrants, three of which have been used for the analyses since the fourth quadrant does not have any data. The details of the quadrants for both the wetting phase and non-wetting phase are listed in
Table 3. The
S and
X values at which these quadrants have been divided are based on having an approximately equal amount of data in them.
Figure 6 and
Figure 7 show the sensitivities of
S,
X, and
A in the three quadrants for the wetting and non-wetting phases, respectively. For the wetting phase (
Figure 6), we observe that saturation dominates all quadrants over connectivity and specific area. In quadrant 2, saturation is significantly more important than both
X and
A. The high importance is likely owing to the residual saturations in this quadrant of the
S-X space and thus at these residual zones, phase saturation becomes the most important parameter. We also observe that in quadrant 3, with high connectivity and high saturation, the relative importance of specific area increases compared to quadrants 1 and 2.
For the non-wetting phase (
Figure 7), connectivity is significantly more important than saturation and specific area in quadrant 1. In quadrants 2 and 3, connectivity and saturation are more important than A. The high importance of connectivity in quadrant 1 can again be due to the residual saturations that lie in this quadrant of the
S-X space. At these residual points, the non-wetting phase becomes trapped, leading to the large magnitude of the phase connectivity in quadrant 1. Like the wetting phase, for the non-wetting phase in quadrant 3, with high saturation and high connectivity, the relative importance of specific area increases as compared to quadrant 1 and 2.
Although specific area is not as important as connectivity or saturation from the above sensitivity results, it still plays an important role in the functional form of relative permeability. Such a qualitative sensitivity analysis of the out-of-bag method is one of its drawbacks.
4. Conclusions
We calculated the specific interfacial area (A) for the secondary drainage cycles and observed that it follows a similar trend as primary drainage and imbibition cycles. The specific interfacial area for the secondary drainage starts at the termination point of the imbibition and then decreases with the increasing saturation of the invading phase due to the reconnection of the snapped of phase in the model.
Relative permeability is most sensitive to saturation for the wetting phase and to phase connectivity for the non-wetting phase. Although specific interfacial area is not as important as the other two predictors (S and X) for relative permeability, it can still aid in the formulation of an accurate state function for relative permeability. The study gives insights into the dependence of relative permeability on the three state parameters S, X, and A using a random forest ensemble to estimate the importance of these parameters. The importance of the parameters studied is qualitative and a more quantitative approach may be useful.