Regular Article
A modified celltocell simulation model to determine the minimum miscibility pressure in tight/shale formations
School of Mining and Petroleum Engineering, Faculty of Engineering, University of Alberta, T6G 1H9 Edmonton, AB, Canada
^{*} Corresponding author: huazhou@ualberta.ca
Received:
14
March
2021
Accepted:
20
May
2021
A new oil–gas Minimum Miscibility Pressure (MMP) calculation algorithm is developed in this work based on the classic celltocell simulation model. The proposed algorithm couples the effects of capillary pressure and confinement in the original celltocell simulation model to predict the oil–gas MMPs in a confined space. Given that the original celltocell algorithm relies on the volume predictions of the reservoir fluids in each cell, a volumetranslated PengRobinson Equation of State (PREOS) is applied in this work for improved accuracy on volume calculations of the reservoir fluids. The robustness of the proposed algorithm is examined by performing the confined MMP calculations for four oil–gas systems. The tieline length extrapolation method is used to determine the oil–gas MMP in confined space. The oil recovery factor calculated by the proposed MMP calculation algorithm is then used to validate the results. First, to achieve stable modeling results for all four examples, a total cell number of 500 is determined by examining the variations in the oil recovery as a function of cell number. Then, by calculating the oil recovery factor near the MMP region, it is found that the MMP determined by tieline length method is slightly lower than the inflection point of the oil recovery curve. Through the case studies, the effects of temperature, pore radius, and injection gas impurity on the confined oil–gas MMP calculations are studied in detail. It is found that the oil–gas MMP is reduced in confined space and the degree of this reduction depends on the pore radius. For all the tested pore radii, the confined MMP first increases and then decreases with an increasing temperature. Furthermore, compared to pure carbon dioxide (CO_{2}) injection, the addition of methane (CH_{4}) in the injection gas increases the oil–gas MMP in confined nanopores. Therefore, it is recommended to control the content of CH_{4} in the injection gas in order to achieve a more efficient gas injection design.
© H. Sun & H. Li, published by IFP Energies nouvelles, 2021
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
The gas injection process has long been applied in the petroleum industry as an effective method of enhancing oil recovery [1]. It offers the advantage of a theoretical 100% local oil recovery efficiency over other enhanced oil recovery methods [2]. More recently, it has also been widely applied in unconventional light to medium oil reservoirs and has led to major oil recovery increases [3, 4]. The oil–gas Minimum Miscibility Pressure (MMP), being the most crucial parameter in a gas injection project, needs to be accurately predicted to maximize the efficiency of a gas injection process. The conventional methods of determining the oil–gas MMP are empirical correlations [5–12], experimental methods [13–18], and numerical simulations [19–22]. The empirical correlations are simple and fast to use, but they are developed with limited experimental data and they could become invalid under certain circumstances. The experimental methods provide reliable and accurate predictions on the oil–gas MMPs under reservoir conditions, but they are very expensive and timeconsuming. Provided that all the empirical correlations and the experimental methods come with inherent drawbacks, numerical simulations have become the most popular methods to offer oil–gas MMP estimations in both conventional and unconventional reservoirs with decent accuracy.
The numerical simulation models currently applied in studying the oil–gas MMP in conventional reservoirs include the one cell simulation algorithm [23], the vanishing InterFacial Tension (IFT) calculation method [17, 24], the original celltocell simulation method [19], and the MultipleMixingCell (MMC) method [22]. One cell simulation model is the simplest model to estimate the oil–gas MMP [23]. In this model, only a single cell is considered, and only forward or backward contacts are performed to determine the oil–gas MMP corresponding to a vaporizing or a condensing drive miscibility. This method is useful when the miscibility is controlled by the gas tieline or the oil tieline [23]. It becomes invalid when the miscibility is controlled by a crossover tieline because the one cell simulation tends to overestimate the oil–gas MMP when the miscibility is induced by a vaporizing/condensing drive mechanism [25]. The vanishing IFT method was originally proposed as a criterion in measuring the oil–gas MMP in experiments [17]. The essence of this method is that the oil–gas MMP is determined when the vapor–liquid IFT reaches zero. Later, researchers modified this method and used the Parachor model [24] combined with a cubic equation of states to calculate the vapor–liquid IFT to determine the theoretical oil–gas MMP [26–29]. The original celltocell simulation method was first proposed by Metcalfe et al. [19]. This method mimics the physical 1D gas displacement process and provides the results of final oil recovery, as well as the oil–gas MMP. This model is slower than the previous two methods (i.e., the one cell simulation and the vanishing IFT model) but offers more reliable results. The MMP determination criteria for this method include tracking the tieline path in a ternary diagram [19], tracking the final oil recovery at different pressures until the oil recovery reaches 97% [30], and tracking the minimum key tieline length at different pressures until it reaches zero [31, 32]. The MMC model proposed by Ahmadi and Johns [22] is a simplified and accelerated model of the original celltocell model [19]. The minimum key tieline length is recorded at each pressure and the last few minimum key tieline lengths are extrapolated to zero, at which the oil–gas MMP is determined [22]. However, in this method, the volumes of each cell and mixing fluids are not calculated. As a result, oil recovery calculations are not made available in this model.
The aforementioned MMP determination methods are all proposed to predict the oil–gas MMP in conventional oil reservoirs. They tend to lose accuracy in calculating the oil–gas MMP in unconventional reservoirs due to the strong capillarity and confinement effects introduced by the nanopores prevalent in unconventional formations. One of the first attempts to reflect the effect of nanopores on the oil–gas MMP prediction is performed by Teklu et al. [33]. They coupled the capillary pressure and critical point shift model [34] in the MMC method to calculate the oil–gas MMP for three oil–gas systems. They found that the oil–gas MMPs are greatly reduced in nanopores compared to the bulk MMPs. It was also concluded that most of the MMP reductions in nanopores are caused by the effect of critical point shift [33]. Other efforts have tried to modify the vanishing IFT calculation model by coupling the effect of nanopores. Teklu et al. [26] coupled the capillary pressure and critical point shift [34] in vapor–liquid equilibrium calculations and calculated vapor–liquid IFT using Parachor model [24]. They found that using the proposed model, the oil–gas MMP decreases with a decreasing pore radius. Similar results have also been observed in the study of Wang et al. [27]. They employed the PerturbedChain Statistical Associating Fluid Theory (PCSAFT) [35] to calculate the vapor–liquid IFT with increasing pressure. They also considered IFT reduction due to the confinement effect. A significant MMP reduction in confined nanopores is observed in their research. Zhang et al. [28] coupled the effect of nanopores in PengRobinson Equation of State (PREOS) [36] and calculated the confined oil–gas MMP using the vanishing IFT model. In addition, the interface thickness is considered in their model by taking account of the twoway mass transfer effects. The proposed model is tested with two crude oil samples and it was found that the oil–gas MMP is reduced in nanopores for both oil samples [28]. Zhang et al. [29] applied a similar methodology and observed a linear increase of the confined oil–gas MMP with an increasing temperature.
The MMC method [22] modified for tight/shale formations is incapable of calculating the oil recovery factors during the gas displacement process. In addition, volume translation is absent in most of PREOSbased confined MMP calculation models. This work aims to develop a comprehensive confined MMP calculation algorithm to predict the oil–gas MMP in confined nanopores. The proposed MMP calculation algorithm is based on the original celltocell simulation model proposed by Metcalfe et al. [19] and the minimum tieline length criterion is applied to determine the confined oil–gas MMP. The celltocell model is selected as the very foundation of this study because this model offers the results of the final oil recovery and the calculated recovery factors can be used to determine the required total cell number in a robust manner to avoid numerical dispersion issues. Also, the fluidmoving mechanism of the celltocell model is more realistic than that of the MMC method. The vapor–liquid equilibrium calculations in confined nanopores are conducted using a volume translated PREOS [36, 37] coupled with capillary pressure and a stateofart critical point shift model [38]. The robustness of the proposed algorithm is examined using four examples. The effects of temperature, pore radius, and injection gas impurity on the confined oil–gas MMP calculations are investigated in detail.
2 Methodology
The mathematical formulations of the thermodynamic model used in this work is presented in this section. A volume translated PREOS [36, 37] is employed in the twophase equilibrium calculations. The capillary pressure and a critical point shift model are coupled with PREOS to reflect the effects of capillarity and confinement yielded by nanopores.
The twophase equilibrium is reached when the fugacityequality condition is satisfied at given pressure, temperature, and mixture compositions. The fugacityequality condition in confined nanopores is given in equation (1) [39],
(1)where f_{ix} and f_{iy} represent the fugacities of component i in the liquid phase and the vapor phase, respectively. P_{l} is the liquidphase pressure, P_{v} is the vaporphase pressure. T is the reservoir temperature, x_{i} is the composition of the liquid phase, y_{i} is the composition of the vapor phase, and n_{c} is the number of components. The component fugacities can be calculated by equation (2) [39],
(2)where P is the phase pressure, ϕ_{i} is the fugacity coefficient of the ith component in any phase, Z is the phase compressibility factor, and A and B are EOS constants. Because of the existence of large capillary pressure, the vaporphase pressure and the liquidphase pressure cannot be assumed equal and can be related using equation (3),
(3)where P_{c} is the capillary pressure of the adjacent liquid and vapor phases. The vapor–liquid capillary pressure can be calculated using the YoungLaplace equation [40] with assumptions of zero contact angle and equal principal radii. The YoungLaplace equation is expressed in equation (4) [40],
(4)where r_{p} is the pore radius, σ represents the IFT between the liquid phase and the vapor phase. The Parachor model [24] is employed to calculate the vapor–liquid IFT given in equation (5),
(5)where P_{chi} is the Parachor constant of component i, ρ_{l} and ρ_{g} are the molar density of the liquidphase and the vaporphase, respectively, and can be calculated by the volume translated PREOS [36, 37].
In addition to capillary pressure, the strong wallmolecule interactions present in nanopores also alter the critical point of the constituting components in reservoir fluids [41]. An analytical solution to express the critical point shift due to the confinement effect was originally proposed by Zarragoicoechea and Kuz [34]. More recently, Tan et al. [38] developed a new critical point shift model based on the one developed by Zarragoicoechea and Kuz [34]. The model proposed by Tan et al. [38] is validated by experiments, making it more representative in describing the critical point shift in confined nanopores. The critical point shift model is given in equation (6) [38]:
(6)where T_{cb} and T_{cp} are the critical temperatures in the bulk conditions and nanopores, respectively, P_{cb} and P_{cp} are the critical pressure in the bulk conditions and nanopores, respectively, and σ_{LJ} is the LennardJones size parameter. It is worth noting that although the capillary pressure and a critical shift model [38] are coupled in the proposed algorithm, this work does not reflect the effect of interactions at the molecular level, such as adsorption, on the confined MMP calculations.
The entire twophase equilibrium calculation algorithm is started with the initialization of the phase equilibrium ratio (k_{i}) using equation (7) [39]:
(7)where k_{i} is the phase equilibrium ratio of the ith component, P_{ci} is the critical pressure of the ith component, P is the system pressure, ω_{i} is the acentric factor of the ith component, T_{ci} is the critical temperature of the ith component, and T is the system temperature. The k_{i} is then used to solve the RachfordRice equation for the phase mole fraction and the compositions of the vapor and the liquid phases. The RachfordRice equation is given in equation (8) [42],
(8)where β_{y} is the vapor phase mole fraction, and z_{i} is the mixture feed composition. The obtained phase mole fraction and the phase compositions are applied to examine the fugacity equality condition with an error tolerance of 10^{−10}. If the fugacity equality condition is reached, then the algorithm stops, and the phase mole fraction and the phase compositions are obtained as outputs. If the fugacity equality condition is not reached, the phase equilibrium ratio is updated using equation (9) and the above procedure is repeated until the fugacity equality condition is reached,
(9)where n represents the current iteration and n + 1 represents the next iteration.
3 Proposed MMP calculation algorithm
The proposed MMP calculation algorithm is based on the original celltocell simulation model developed by Metcalfe et al. [19]. A detailed stepbystep procedure of performing the confined MMP calculations using the proposed algorithm is shown as follows:

The reservoir temperature, compositions of injection gas and reservoir oil, cell volume (v_{c}), Gas/Oil Ratio (GOR), and total cell number (N_{c}) are specified to initiate the algorithm. It is worth noting that the cell volume is not a significant parameter and is assumed as 1 cm^{3}. The GOR is set to be 0.3. It is shown by Zhao et al. [31] that a larger GOR may result in poor accuracy of oil recovery calculations, while a smaller GOR leads to a significant increase of computational cost. Then the shifted critical properties of each constituting component are calculated using the critical point shift model [38].

The total volume of injection is set to be 1.2 Pore Volume (PV) as suggested by Metcalfe et al. [19]. This is also a standard volume of injection in gas and oil industry for a gas injection design. This volume of injection ensures an optimum efficiency of a gas injection [19]. Then the total batch number (N_{b}) can be calculated using equation (10),
(10)where N_{b} is the total batch number, N_{c} is the total cell number, and GOR is the Gas/Oil Ratio.

The volume of the first batch of injection gas is calculated as GOR × v_{c}. Then the first contact is made by mixing the first batch of injection gas with reservoir oil. It is assumed that all the cells are initially filled with reservoir oil and perfect mixing is reached with each contact. Therefore, the feed compositions can be calculated using the volume of the injection gas and the reservoir oil. Then the compositions and the phase volumes of the resulting mixture can be calculated using the twophase equilibrium calculation algorithm coupled with the capillary pressure and critical point shift model [38]. The phase volumes are predicted using the volume translated PREOS [36, 37] for better accuracy. In addition, the tieline length of this contact is calculated using equation (11) and recorded [22],
(11) where TL is the TieLine length, x_{i} and y_{i} represent the mole fractions of component i in the equilibrium liquid phase and the equilibrium vapor phase, respectively.

The moving strategy used in the proposed algorithm is to move the excess fluid. If the resulting mixture after the twophase equilibrium calculation is singlephase, then the excess gas or oil is moved to the next cell. If the resulting mixture is a twophase system, then the excess gas is moved to the next cell first. If after all the gas is moved to the next cell and the volume of the remaining oil is still larger than the cell volume, then the excess oil is also moved to the next cell [19].

The overall feed compositions of the fluids in cell 2 are calculated and the twophase equilibrium calculation is conducted using the calculated overall feed compositions. The tieline length is also recorded for this contact. The excess fluid is then moved to cell 3. This step is repeated until the last cell.

After one batch of gas injection is completed, another batch of gas is injected into cell 1 until all the 1.2 PV of gas is injected. The fluid moved from the last cell is treated as the final oil production. The final oil recovery factor is calculated using equation (12) [31],
(12)where RF_{1.2} is the oil recovery factor after 1.2 PV gas injection, V_{r} is the volume of the oil recovered from the last cell at 1 atm and 298 K, V_{o} is the volume of the original oil in place at 1 atm and 298 K. The minimum tieline length among the tieline lengths of all contacts is recorded.

The above steps are repeated at different pressures and the minimum tieline length is recorded for each pressure. The minimum tieline length is reduced at a higher pressure and the MMP can be theoretically found at the pressure where the minimum tieline length reaches zero. However, it is difficult to calculate the tieline length when the pressure approaches the MMP [33]. This issue can be effectively addressed by extrapolating the minimum tieline length at the last few pressures [22]. A powerlaw extrapolation as shown in equation (13) is applied [22],
(13)where n is the exponential parameter, a is the slope, and b is the intercept of the yaxis. It is recommended that the square of the correlation coefficient should exceed 0.999 when finalizing the parameters n, a, and b [22].

The confined oil–gas MMP determined by the tieline length method is then validated by the oil recovery factors calculated by the proposed algorithm. At a constant temperature, the oil recovery factor at 1.2 PV gas injection is calculated at different pressures. The MMP determined by the tieline length method is considered reasonable if it is near the inflection point of the oil recovery curve.
Compared to the original celltocell simulation model [19] which can only estimate the oil–gas MMP under bulk conditions, the proposed algorithm is capable of predicting the oil–gas MMP in a confined space. Also, the proposed algorithm employs a volume translated PREOS [36, 37] to improve the accuracy of phase volume predictions. Compared to the MMC method proposed by Ahmadi and Johns [22], the celltocell simulation model [19] on which this work is based has the advantage of calculating the final oil recovery at each pressure which can be used to determine the required total cell number in order to achieve stable MMP results. The calculated oil recovery factor can also be used to validate the MMP results determined by the tieline length method. Therefore, the calculated final oil recovery profile is used as a secondary validation of the proposed algorithm. This can effectively improve the reliability of the confined MMP calculations and is not featured in the previous confined MMP calculation algorithms. A schematic of the celltocell simulation model and an overall procedure of conducting the confined MMP calculations using the proposed algorithm can be visualized in the flowchart shown in Figure 1.
Fig. 1 (a) A schematic of the celltocell simulation model and (b) a workflow of conducting the confined oil–gas MMP calculations using the proposed algorithm. 
4 Results and discussions
This section presents the confined oil–gas MMP calculations for four oil–gas systems using the proposed algorithm. The effects of temperature, pore radius, and injection gas impurity on the confined MMP calculations are also studied.
4.1 Summary of the tested oil–gas systems
A total of four oil–gas systems are tested in this work. The first three fluid systems are simple synthetic oil samples being displaced by CO_{2}, CH_{4,} and their mixture. The detailed compositions of the injection gas and the oil for the first three fluid systems are summarized in Table 1.
The physical properties of the gas/oil constituting components for the first three oil–gas systems and their Binary Interactive Parameters (BIPs) with CO_{2} used in the EOS model are given in Table 2.
The BIPs between hydrocarbon components are assumed to be zeros as they are negligible compared to hydrocarbon–CO_{2} BIPs. The last fluid system is a real dead crude oil sample being displaced by pure CO_{2} [44]. The BIPs of the hydrocarbon components in this oil sample are tuned by matching the calculated confined MMP with the experimental data and the detailed procedure can be found in the work of Sun and Li [45]. The physical properties of the constituting components and their BIPs with CO_{2} for the fourth oil–gas system are given in Table 3.
The confined oil–gas MMPs of the above four fluid systems at different temperatures and pore radii are calculated using the proposed algorithm. The obtained confined MMPs are then analyzed to study the effects of temperature, pore radius, and injection gas impurity on the confined MMP calculations.
4.2 Determining the total cell number
The total cell number and GOR are the two factors that might cause severe numerical dispersion issues. The GOR in this work is set to be 0.3 as suggested by Metcalfe et al. [19]. The total cell number, however, has not received a widely accepted criterion that can eliminate the numerical dispersion. In this work, the final oil recovery is calculated with different total cell numbers until the final oil recovery flattens out with an increasing total cell number. The appropriate total cell number can be selected to be the one at which the stable oil recovery is maintained in both bulk conditions and in confined nanopores. The final oil recoveries as a function of the total cell number for all four oil–gas systems are shown in Figure 2.
Fig. 2 The effect of the total cell number on the final oil recovery under bulk conditions and in confined nanopores for the (a) oil–gas system 1, (b) oil–gas system 2, (c) oil–gas system 3, and (d) oil–gas system 4. 
It can be seen from Figure 2 that the calculated final oil recovery increases with an increasing total cell number in both bulk and confined spaces. The final oil recovery for the first three fluids systems becomes stable at a total cell number of 150. The fourth fluid system, however, requires 450 to reach a stable final oil recovery. This can be attributed to the presence of the heavy components in the crude oil sample of the fourth fluid system. The injected CO_{2} needs to be enriched with additional contacts before it is able to reach miscibility with the crude oil. In order to guarantee reliable and accurate results, this work applies a total cell number of 500 in all the confined oil–gas MMP calculations. Figure 3 shows the change of the resulting gas and oil compositions after each contact at different cell numbers for the fluid system 3 as an example to visualize the gas enrichment process through multiple contacts in bulk conditions and in a confined space.
Fig. 3 The mole fractions of the constituting components in oil phase and gas phase for the oil–gas system 3 at a temperature of 306.15 K and a pressure of (a) 220 bar in bulk conditions, (b) 220 bar in a confined space, (c) 250 bar in bulk conditions, and (d) 250 in a confined space. 
The solid lines and the dashed lines in Figure 3 represent the compositions of the liquid phase and the gas phase, respectively. Along the displacing path, the injection gas is enriched with C_{4} and C_{10}. It is seen that the gas enrichment path in bulk conditions is different from that in a confined space. In general, the gas is enriched with more hydrocarbon components in a confined space than in bulk conditions at a fixed pressure. For the gas enrichment process in a confined space, the calculated MMP at 306.15 K with a pore radius of 10 nm is 251.02 bar. From Figure 3, it is found that with an increasing pressure, the injection gas is enriched with more hydrocarbon components. Also, when the pressure approaches the MMP, the oil and gas compositions become almost identical. This also proves the correctness of the proposed algorithm.
4.3 Further validation of MMP using the calculated oil recovery factors
As described before, the calculated oil recovery from the proposed MMP calculation algorithm is used to validate the confined oil–gas MMP determined by the tieline length method. As per the slimtube test protocols, the MMP can be determined as the inflection point of the oil recovery curve [14]. This criterion is also used in this work to validate the MMP results obtained from the tieline length method. It is worth noting that under most circumstances, the Wilson equation [46] (e.g., Eq. (7)) is applied to provide a reliable initialization of the phase equilibrium ratios (k_{i}). In the nearMMP region or the region above MMP, other k_{i} initialization strategies are required including k_{i}/3 and [39]. Figure 4 shows the oil recovery curves calculated for all four fluid systems at different temperatures with a pore radius of 2 nm. In Figure 4, the confined MMPs determined using the tieline length method are also highlighted using vertical lines.
Fig. 4 Calculated oil recovery as a function of pressure for the (a) fluid system 1, (b) fluid system 2, (c) fluid system 3, and (d) fluid system 4 at different temperatures with a pore radius of 2 nm. 
It can be seen from Figure 4 that the oil recovery factor increases with an increasing pressure. Above a certain pressure, the increase rate of oil recovery is reduced dramatically. The inflection point can be determined as the oil–gas MMP [14]. It is found from Figure 4 that for all the four fluid systems, the confined MMPs determined by the tieline length method are slightly lower than the inflection points of the oil recovery curve for all tested temperatures. Nevertheless, the MMPs determined by the tieline length method are very close to the inflection point of the oil recovery curve. We also compare the confined MMPs calculated using the proposed algorithm with the ones calculated using the confined MMC method [45]. Figure 5 shows the comparison results in a parity chart. It can be seen from Figure 5 that the confined MMPs calculated by the proposed algorithm and the MMC method [45] are in an excellent agreement. This is expected because both methods employ the tieline length extrapolation method to determine the confined MMPs. Furthermore, for this oil sample, the experimentally measured CO_{2}–oil MMP in the tested shale cores is 123 bar at 326.15 K [44]. The confined MMP calculated using the proposed algorithm is 122.38 bar under experimental conditions. This demonstrates that the proposed algorithm is able to reproduce the experimental results using the properly tuned BIPs. The advantage of the proposed algorithm over the previous one is that the proposed algorithm offers the opportunity to validate the determined MMPs using the oil recovery factors. This feature of the secondary validation is crucial because it allows this simulation study to mimic the realistic gas displacement process in a slimtube experiment and hence significantly improves the reliability of the obtained simulation results.
Fig. 5 Confined MMPs calculated by the proposed algorithm and the previously proposed MMC algorithm [45]. 
4.4 Effects of temperature and pore radius on the confined MMP calculations
It is observed in many studies that the oil–gas MMP is heavily dependent on the reservoir temperature and pore size [26, 33, 45]. This work calculates the oil–gas MMP in confined nanopores with different pore sizes at a wide range of temperatures to investigate the effects of temperature and pore radius on confined MMP calculations. The results are incredibly similar to our previous study [45]. The calculated oil–gas MMP in confined nanopores first increases and then decreases with an increasing temperature. It is intriguing to notice that this trend has also been introduced for bulk oil–gas MMP [47]. The reason for the MMP decreases after a certain temperature is that at an elevated temperature, the liquid phase in the reservoir becomes more vaporphase like. Hence, the oil–gas MMP becomes the miscibility pressure between one vapor phase and one vaporphaselike liquid phase. Moreover, a smaller pore radius clearly results in a decreasing oil–gas MMP. The MMP decreases in the confined nanopores from the bulk MMP significantly when the pore radius is less than 5 nm, while it becomes almost constant when the pore radius is larger than 20 nm. The above findings are in high consistency with the ones obtained in Sun and Li [45] which applied a modified MMC method in calculating the confined MMPs.
4.5 Injection gas impurity effect on the confined MMP calculations
Another advantage of the proposed MMP algorithm is that it can be used to conveniently investigate the effect of injection gas compositions on the gas–oil MMP. As an example, here we employ the proposed algorithm to examine the effect of CH_{4} contamination in CO_{2} on the gas–oil MMP. Due to the difficulty of acquiring pure CO_{2} for the injection, it is a common practice in the field to reinject the produced CH_{4} together with CO_{2} back into the reservoir for EOR purposes [29, 48]. The previous calculation results for the first and the second fluid systems clearly show that the confined MMP is generally increased when injecting a CH_{4}–CO_{2} mixture as compared to pure CO_{2} injection. Herein, we calculate the confined MMP of the first and the second fluids systems being displaced by CO_{2}–CH_{4} mixture with various CH_{4} content. Figure 6 shows the calculation results.
Fig. 6 Calculated confined oil–gas MMP as a function of the CH_{4} mole fraction in the injection gas. 
The calculation results in Figure 6 suggest that the calculated oil–gas MMP in confined nanopores increases with an increasing CH_{4} content. The MMP increasing trend appears to be almost linear, which is similar to the MMP predictions in bulk conditions [49]. Such observation seems to be true for all the tested pore radii. Therefore, in a real gas injection practice in tight/shale reservoirs, the content of CH_{4} in the injection gas should be limited in order to achieve a lower oil–gas MMP. Through this example, it is demonstrated that the proposed algorithm can be readily used to study the effect of injection gas impurity on the confined MMP calculations.
5 Conclusion
We develop a new oil–gas MMP calculation algorithm based on the original celltocell simulation model [19] to predict the oil–gas MMP in confined nanopores with decent reliability and accuracy. The results of the final oil recovery are utilized to determine the required total cell number to avoid numerical dispersion and to validate the confined MMPs determined by the tieline length method. The capillary pressure and a critical point shift model [38] are coupled with the conventional twophase equilibrium calculation algorithm to reflect the capillarity and confinement effect in nanopores. A volume translated PREOS [36, 37] is employed to provide accurate phase densities calculations. By conducting confined oil–gas MMP calculations for four oil–gas systems at different temperatures and pore radii, we can reach the following conclusions:
A total cell number of 500 is determined as the cell number that can give stable simulation results. Compared to a light to medium oil sample, the oil sample containing more heavy components requires a larger total cell number to reach a stable oilrecoveryfactor calculation result. This is because the injection gas requires more contacts to be enriched when displacing a heavier oil as compared to the displacement of a lighter oil.
The confined MMP determined by the tieline length method can be further validated by the calculated oil recovery curves. It is found that the confined MMPs determined by the tieline length method are slightly lower than the inflection point of the oil recovery curve. The calculated confined MMP corresponds to the pressure at which the increase of the oil recovery factor starts to slow down.
The effects of the temperature and pore radius on the confined gas–oil MMP calculations discussed in this work are in an excellent agreement with our previous study [45] where a modified MMC method was applied. The confined oil–gas MMP first increases and then decreases with an increasing temperature. This is due to the phase behavior limit where the oil phase becomes a vaporlike phase or vapor phase at an elevated temperature. Also, the oil–gas MMP in confined nanopores is reduced. The reduction of the confined MMP becomes significant when the pore radius is less than 10 nm. When the pore radius is larger than 20 nm, the oil–gas MMP becomes almost constant and is close to the bulk MMP.
The injection gas impurity has a strong effect on the confined oil–gas MMP. The simulation results using the proposed MMC algorithm suggest that the confined MMP increases with an increasing CH_{4} content in the injection gas.
Acknowledgments
The authors acknowledge a Discovery Grant from the Natural Sciences and Engineering Research Council of Canada (NSERC) to H. Li (Grant No.: NSERC RGPIN202004571).
References
 Belhaj H., Abukhalifeh H., Javid K. (2013) Miscible oil recovery utilizing N_{2} and/or HC gases in CO_{2} injection, J. Pet. Sci. Eng. 111, 144–152. [Google Scholar]
 Elsharkawy A.M., Poettmann F.H., Christiansen R.L. (1996) Measuring CO_{2} minimum miscibility pressures: Slimtune or risingbubble method? Energy Fuels 10, 443–449. [Google Scholar]
 Yu W., Lashgari H., Wu K., Sepehrnoori K. (2015) CO_{2} injection for enhanced oil recovery in bakken tight oil reservoirs, Fuel 159, 354–363. [Google Scholar]
 Lashgari H., Sun A., Zhang T., Pope G., Lake L. (2019) Evaluation of carbon dioxide storage and miscible gas EOR in shale oil reservoirs, Fuel 241, 1223–1235. [Google Scholar]
 Holm L., Josendal V. (1974) Mechanisms of oil displacement by carbon dioxide, J. Pet. Technol. 26, 12, 1427–1438. [Google Scholar]
 Lee J.I. (1979) Effectiveness of carbon dioxide displacement under miscible and immiscible conditions, Report RR40, Petroleum Recovery Institute, Calgary. [Google Scholar]
 Mungan N. (1981) Carbon dioxide floodingfundamentals, J. Can. Pet. Technol. 20, 01, 87–92. [Google Scholar]
 Orr F., Jensen C. (1984) Interpretation of pressurecomposition phase diagrams for CO_{2}/crudeoil systems, Soc. Pet. Eng. J. 24, 05, 485–497. [Google Scholar]
 Zhang T., Li Y., Sun S. (2019) Phase equilibrium calculations in shale gas reservoirs, Capillarity 2, 1, 8–16. [Google Scholar]
 ZareNezhad B. (2016) A new correlation for predicting the minimum miscibility pressure regarding the enhanced oil recovery processes in the petroleum industry, Pet. Sci. Technol. 34, 1, 56–62. [Google Scholar]
 Li H., Qin J., Yang D. (2012) An improved CO_{2}–oil minimum miscibility pressure correlation for live and dead crude oils, Ind. Eng. Chem. Res. 51, 8, 3516–3523. [Google Scholar]
 Choubineh A., Helalizadeh A., Wood D.A. (2019) Estimation of minimum miscibility pressure of varied gas compositions and reservoir crude oil over a wide range of conditions using an artificial neural network model, Adv. GeoEnergy Res. 3, 1, 52–66. [Google Scholar]
 Rathmell J.J., Stalkup F.I., Hassinger R.C. (1971) A laboratory investigation of miscible displacement by carbon dioxide, in: Fall Meeting of the Society of Petroleum Engineers of AIME, New Orleans, Louisiana, USA, October 3–6, 1971, Society of Petroleum Engineers. [Google Scholar]
 Yellig W. (1982) Carbon dioxide displacement of a West Texas reservoir oil, SPE J. 22, 06, 805–815. [Google Scholar]
 Christiansen R.L., Kim H. (1986) Apparatus and method for determining the minimum miscibility pressure of a gas in a liquid, Energy Fuels 10, 443–449. [Google Scholar]
 Randall T., Bennion D. (1988) Recent developments in slim tube testing for HydroCarbonMiscible Flood (HCMF) solvent design, J. Can. Pet. Technol. 27, 06, 33–44. [Google Scholar]
 Rao D. (1997) A new technique of vanishing interfacial tension for miscibility determination, Fluid Phase Equilib. 139, 1–2, 311–324. [Google Scholar]
 Nguyen P., Mohaddes D., Riordon J., Fadaei H., Lele P., Sinton D. (2015) Fast fluorescencebased microfluidic method for measuring minimum miscibility pressure of CO_{2} in crude oils, Anal. Chem. 87, 6, 3160–3164. [Google Scholar]
 Metcalfe R.S., Fussell D.D., Shelton J.L. (1973) A multicell equilibrium separation model for the study of multiple contact miscibility in richgas drives, SPE J. 13, 3, 147–155. [Google Scholar]
 Johns R., Orr F. (1996) Miscible gas displacement of multicomponent oils, SPE J. 1, 01, 39–50. [Google Scholar]
 Wang Y., Orr F. (1997) Analytical calculation of minimum miscibility pressure, Fluid Phase Equilib. 139, 1–2, 101–124. [Google Scholar]
 Ahmadi K., Johns R. (2011) Multiplemixingcell method for MMP calculations, SPE J. 16, 04, 733–742. [Google Scholar]
 Neau E., Avaullée L., Jaubert J.N. (1996) A new algorithm for enhanced oil recovery calculations, Fluid Phase Equilib. 117, 265–272. [Google Scholar]
 Weinaug C., Katz D. (1943) Surface tensions of methanepropane mixtures, Ind. Eng. Chem. 35, 2, 239–246. [Google Scholar]
 Jaubert J.N., Wolff L., Neau E., Avaullée L. (1998) A very simple multiple mixing cell calculation to compute the minimum miscibility pressure whatever the displacement mechanism, Ind. Eng. Chem. Res. 37, 4854–4859. [Google Scholar]
 Teklu T., Alharthy N., Kazemi H., Yin X., Graves R., AlSumaiti A. (2014) Phase behavior and minimum miscibility pressure in nanopores, SPE Reserv. Eval. Eng. 17, 03, 396–403. [Google Scholar]
 Wang S., Ma M., Chen S. (2016) Application of PCSAFT equation of state for CO_{2} minimum miscibility pressure prediction in nanopores, in: SPE Improved Oil Recovery Conference, Tulsa, Oklahoma, USA, April 11–13, 2016, Society of Petroleum Engineers. [Google Scholar]
 Zhang K., Jia N., Zeng F., Luo P. (2017) A new diminishing interface method for determining the minimum miscibility pressures of light oil–CO_{2} systems in bulk phase and nanopores, Energy Fuels 31, 11, 12021–12034. [Google Scholar]
 Zhang K., Jia N., Li S. (2017) Exploring the effects of four important factors on oil–CO_{2} interfacial properties and miscibility in nanopores, RSC Adv. 7, 85, 54164–54177. [Google Scholar]
 Jaubert J.N., Arras L., Neau E., Avaullée L. (1998) Properly defining the classical vaporizing the condensing mechanisms when a gas is injected into a crude oil, Ind. Eng. Chem. Res. 37, 4860–4869. [Google Scholar]
 Zhao G., Adidharma H., Towler B., Radosz M. (2006) Using a multiplemixingcell model to study minimum miscibility pressure controlled by thermodynamic equilibrium tie lines, Ind. Eng. Chem. Res. 45, 23, 7913–7923. [Google Scholar]
 Yang F., Yu P., Zhang X. (2020) Multiplemixingcell model for calculation of minimum miscibility pressure controlled by tieline length, Geofluids, 2020, 1–8. [Google Scholar]
 Teklu T.W., Alharthy N., Kazemi H., Yin X., Graves R.M. (2014) Vanishing interfacial tension algorithm for MMP determination in unconventional reservoirs, in: SPE Western North American and Rocky Mountain Joint Regional Meeting, Denver, Colorado, USA, April 17–18, 2014, Society of Petroleum Engineers. [Google Scholar]
 Zarragoicoechea G., Kuz V. (2004) Critical shift of a confined fluid in a nanopore, Fluid Phase Equilib. 220, 1, 7–9. [Google Scholar]
 Gross J., Sadowski G. (2002) Application of the perturbedchain SAFT equation of state to associating systems, Ind. Eng. Chem. Res. 41, 22, 5510–5515. [Google Scholar]
 Peng D., Robinson D. (1976) A new twoconstant equation of state, Ind. Eng. Chem. Fundam. 15, 1, 59–64. [Google Scholar]
 Abudour A., Mohammad S., Robinson R., Gasem K. (2013) Volumetranslated PengRobinson equation of state for liquid densities of diverse binary mixtures, Fluid Phase Equilib. 349, 37–55. [Google Scholar]
 Tan S., Qiu X., Dejam M., Adidharma H. (2019) Critical point of fluid confined in nanopores: Experimental detection and measurement, J. Phys. Chem. C 123, 15, 9824–9830. [Google Scholar]
 Whitson C., Brulé M. (2000) Phase behavior, Henry L. Doherty Memorial Fund of AIME, Society of Petroleum Engineers, Richardson, Texas. [Google Scholar]
 Young T. (1805) An essay on the cohesion of fluids, Philos. Trans. R. Soc. Lond. 95, 65–87. [Google Scholar]
 Qiu X., Tan S., Dejam M., Adidharma H. (2019) Experimental study on the criticality of a methane/ethane mixture confined in nanoporous media, Langmuir 35, 36, 11635–11642. [Google Scholar]
 Rachford H., Rice J. (1952) Procedure for use of electronic digital computers in calculating flash vaporization hydrocarbon equilibrium, J. Pet. Tech. 4, 193. [Google Scholar]
 Quayle O.R. (1953) The parachors of organic compounds. An interpretation and catalogue, Chem. Rev. 53, 3, 439–589. [Google Scholar]
 Zhang K., Gu Y. (2015) Two different technical criteria for determining the Minimum Miscibility Pressures (MMPs) from the slimtube and coreflood tests, Fuel 161, 146–156. [Google Scholar]
 Sun H., Li H. (2021) Minimum miscibility pressure determination in confined nanopores considering pore size distribution of tight/shale formations, Fuel 286, 119450. [Google Scholar]
 Wilson G.M. (1969) A modified RedlichKwong equation of state, application to general physical data calculations, in: 65th National AIChE Meeting, Cleveland, May 4–7, 1969. [Google Scholar]
 Yuan H., Johns R. (2005) Simplified method for calculation of minimum miscibility pressure or enrichment, SPE J. 10, 04, 416–425. [Google Scholar]
 Zick A.A. (1986) A combined condensing/vaporizing mechanism in the displacement of oil by enriched gases, in: SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, USA, October 5–8, 1986, Society of Petroleum Engineers. [Google Scholar]
 Jin L., Pekot L.J., Hawthorne S.B., Gobran B., Greeves A., Bosshart N.W., Jiang T., Hamling J.A., Gorecki C.D. (2016) Impact of CO_{2} impurity on MMP and oil recovery performance of the bell creek oil field, Energy Procedia 114, 6997–7008. [Google Scholar]
All Tables
All Figures
Fig. 1 (a) A schematic of the celltocell simulation model and (b) a workflow of conducting the confined oil–gas MMP calculations using the proposed algorithm. 

In the text 
Fig. 2 The effect of the total cell number on the final oil recovery under bulk conditions and in confined nanopores for the (a) oil–gas system 1, (b) oil–gas system 2, (c) oil–gas system 3, and (d) oil–gas system 4. 

In the text 
Fig. 3 The mole fractions of the constituting components in oil phase and gas phase for the oil–gas system 3 at a temperature of 306.15 K and a pressure of (a) 220 bar in bulk conditions, (b) 220 bar in a confined space, (c) 250 bar in bulk conditions, and (d) 250 in a confined space. 

In the text 
Fig. 4 Calculated oil recovery as a function of pressure for the (a) fluid system 1, (b) fluid system 2, (c) fluid system 3, and (d) fluid system 4 at different temperatures with a pore radius of 2 nm. 

In the text 
Fig. 5 Confined MMPs calculated by the proposed algorithm and the previously proposed MMC algorithm [45]. 

In the text 
Fig. 6 Calculated confined oil–gas MMP as a function of the CH_{4} mole fraction in the injection gas. 

In the text 