A modified cell-to-cell simulation model to determine the minimum miscibility pressure in tight/shale formations

A new oil–gas MinimumMiscibility Pressure (MMP) calculation algorithm is developed in this work based on the classic cell-to-cell simulation model. The proposed algorithm couples the effects of capillary pressure and confinement in the original cell-to-cell simulation model to predict the oil–gas MMPs in a confined space. Given that the original cell-to-cell algorithm relies on the volume predictions of the reservoir fluids in each cell, a volume-translated Peng-Robinson Equation of State (PR-EOS) 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 tie-line 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 tie-line 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 (CO2) injection, the addition of methane (CH4) in the injection gas increases the oil–gas MMP in confined nanopores. Therefore, it is recommended to control the content of CH4 in the injection gas in order to achieve a more efficient gas injection design.


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 oilgas 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][6][7][8][9][10][11][12], experimental methods [13][14][15][16][17][18], and numerical simulations [19][20][21][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 time-consuming. 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 cell-to-cell simulation method [19], and the Multiple-Mixing-Cell (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 tie-line or the oil tie-line [23]. It becomes invalid when the miscibility is controlled by a crossover tie-line 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][27][28][29]. The original cell-to-cell simulation method was first proposed by Metcalfe et al. [19]. This method mimics the physical 1-D 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 tie-line 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 tie-line 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 cell-to-cell model [19]. The minimum key tie-line length is recorded at each pressure and the last few minimum key tie-line 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 Perturbed-Chain Statistical Associating Fluid Theory (PC-SAFT) [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 Peng-Robinson Equation of State (PR-EOS) [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 two-way 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 PR-EOS-based 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 cellto-cell simulation model proposed by Metcalfe et al. [19] and the minimum tie-line length criterion is applied to determine the confined oil-gas MMP. The cell-to-cell 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 fluid-moving mechanism of the cell-to-cell model is more realistic than that of the MMC method. The vapor-liquid equilibrium calculations in confined nanopores are conducted using a volume translated PR-EOS [36,37] coupled with capillary pressure and a state-of-art 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.

Methodology
The mathematical formulations of the thermodynamic model used in this work is presented in this section. A volume translated PR-EOS [36,37] is employed in the two-phase equilibrium calculations. The capillary pressure and a critical point shift model are coupled with PR-EOS to reflect the effects of capillarity and confinement yielded by nanopores.
The two-phase equilibrium is reached when the fugacity-equality condition is satisfied at given pressure, temperature, and mixture compositions. The fugacity-equality condition in confined nanopores is given in equation (1) [39], 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 liquid-phase pressure, P v is the vapor-phase 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], ln f i 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 vapor-phase pressure and the liquid-phase pressure cannot be assumed equal and can be related using equation (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 Young-Laplace equation [40] with assumptions of zero contact angle and equal principal radii. The Young-Laplace equation is expressed in equation (4) [40], where r p is the pore radius, r 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), where P chi is the Parachor constant of component i, q l and q g are the molar density of the liquid-phase and the vaporphase, respectively, and can be calculated by the volume translated PR-EOS [36,37]. In addition to capillary pressure, the strong wallmolecul 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]: 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 r LJ is the Lennard-Jones 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 two-phase equilibrium calculation algorithm is started with the initialization of the phase equilibrium ratio (k i ) using equation (7) [39]: 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, x 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 Rachford-Rice equation for the phase mole fraction and the compositions of the vapor and the liquid phases. The Rachford-Rice equation is given in equation (8) [42], where b 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, where n represents the current iteration and n + 1 represents the next iteration.

Proposed MMP calculation algorithm
The proposed MMP calculation algorithm is based on the original cell-to-cell simulation model developed by Metcalfe et al. [19]. A detailed step-by-step procedure of performing the confined MMP calculations using the proposed algorithm is shown as follows: 1. 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]. 2. 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), where N b is the total batch number, N c is the total cell number, and GOR is the Gas/Oil Ratio.
3. 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 two-phase equilibrium calculation algorithm coupled with the capillary pressure and critical point shift model [38]. The phase volumes are predicted using the volume translated PR-EOS [36,37] for better accuracy. In addition, the tie-line length of this contact is calculated using equation (11) and recorded [22], where TL is the Tie-Line length, x i and y i represent the mole fractions of component i in the equilibrium liquid phase and the equilibrium vapor phase, respectively.
4. The moving strategy used in the proposed algorithm is to move the excess fluid. If the resulting mixture after the two-phase equilibrium calculation is singlephase, then the excess gas or oil is moved to the next cell. If the resulting mixture is a two-phase 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]. 5. The overall feed compositions of the fluids in cell 2 are calculated and the two-phase equilibrium calculation is conducted using the calculated overall feed compositions. The tie-line length is also recorded for this contact. The excess fluid is then moved to cell 3. This step is repeated until the last cell. 6. 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], 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 tie-line length among the tie-line lengths of all contacts is recorded.
7. The above steps are repeated at different pressures and the minimum tie-line length is recorded for each pressure. The minimum tie-line length is reduced at a higher pressure and the MMP can be theoretically found at the pressure where the minimum tie-line length reaches zero. However, it is difficult to calculate the tie-line length when the pressure approaches the MMP [33]. This issue can be effectively addressed by extrapolating the minimum tie-line length at the last few pressures [22]. A power-law extrapolation as shown in equation (13) is applied [22], where n is the exponential parameter, a is the slope, and b is the intercept of the y-axis. It is recommended that the square of the correlation coefficient should exceed 0.999 when finalizing the parameters n, a, and b [22].
8. The confined oil-gas MMP determined by the tie-line 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 tie-line length method is considered reasonable if it is near the inflection point of the oil recovery curve.
Compared to the original cell-to-cell 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 PR-EOS [36,37] to improve the accuracy of phase volume predictions. Compared to the MMC method proposed by Ahmadi and Johns [22], the cell-to-cell 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 tie-line 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 cell-to-cell 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.

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.

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.

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.
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 Table 1. The compositions of the injection gases and the oils in the first three oil-gas systems [31,32]. Oil-gas system 1 Oil-gas system 2 Oil-gas system 3    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. 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.

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 tie-line length method. As per the slim-tube 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 tie-line 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 near-MMP region or the region above MMP, other k i initialization strategies are required including k i /3 and ffiffiffiffi k i 3 p [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 tie-line length method are also highlighted using vertical lines. 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 tie-line 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 tie-line 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 slim-tube experiment and hence significantly improves the reliability of the obtained simulation results.

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 vapor-phase like. Hence, the oil-gas MMP becomes the miscibility pressure between one vapor phase and one vapor-phase-like 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.

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.
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.

Conclusion
We develop a new oil-gas MMP calculation algorithm based on the original cell-to-cell 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 tie-line length method. The capillary pressure and a critical point shift model [38] are coupled with the conventional two-phase equilibrium calculation algorithm to reflect the capillarity and confinement effect in nanopores. A volume translated PR-EOS [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 oil-recovery-factor 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 tie-line length method can be further validated by the calculated oil recovery curves. It is found that the confined MMPs determined by the tie-line 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 vapor-like 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.