Molecular dynamics simulation to estimate minimum miscibility pressure for oil with pure and impure CO2

CO2 enhanced oil recovery (CO2-EOR) has gained great attention worldwide for increasing oil production and reducing greenhouse emissions. In CO2-EOR project, the minimum miscibility pressure (MMP) is the key factor to determine the possibility of CO2 miscible-phase displacement. In this study, we used molecular dynamics (MD) simulation to estimate the MMP of pure and impure CO2 with oil by obtaining interfacial tensions (IFTs) based on the vanishing interfacial tension technique. To get accurate IFTs, modified Lorentz-Berthelot rules were developed for all simulations. Then, we performed MD simulations for CO2 + n-hexane binary system, CO2 + n-hexane + n-decane ternary system and complex CO2 + crude oil system to calculate the MMPs. Firstly, MMPs predicted by MD simulation show good agreement with previous experimental data. Secondly, the effects of injected gas components, oil components and reservoir temperature on the MMP are investigated. The results demonstrate that methane and N2 will increase the MMP of CO2-oil system significantly; the effects of a certain oil component on the MMP depend on the properties of crude oil, generally speaking, heavy oil components will increase the MMP while light oil components decrease it; temperature will increase the MMP of CO2-oil system, while it will decrease the MMP of N2-oil system, which is due to their different solubility in oil. Thirdly, sensitivity analysis is conducted to evaluate the effects of those influence parameters on the MMP of CO2-oil system, finding that N2 concentration in CO2 and reservoir temperature have more significant impacts on the MMP.


Introduction
Enhanced oil recovery (EOR) has gained considerable attention worldwide, since it can boost the amount of crude oil extracted from the oil reservoir. Among the EOR techniques, CO 2 injection is one of the most promising approaches because it could achieve high oil recovery rate and simultaneously reduce the greenhouse emission by sequestrating CO 2 into depleted oil reservoirs [1][2][3]. To reach the maximum oil recovery, the miscible condition between CO 2 and oil is expected to take place in the reservoir, which requires the operation pressure higher than the minimum miscibility pressure (MMP) [4]. The MMP is defined as the lowest pressure at which the reservoir oil and injected gas obtain 'multiple contact miscibility' at the reservoir temperature [5]. Thus, it is very crucial to know the MMP of CO 2 -oil system before employing CO 2 injection technique. Furthermore, nitrogen and methane are common impurities in CO 2 injection; hence, the MMP for oil with CO 2 , nitrogen, methane and their mixture also are important [6][7][8].
To measure MMP, the slim tube test [9] and the rising bubble apparatus test [10] are the most commonly used experimental methods. Specifically, the slim tube test, involving fluid flow in the sand-packed slim tube, Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. has become the industry standard to determine MMP for decades. However, it is slow and expensive, and lacks a fundamental physiochemical definition of MMP [5,11,12]. The rising bubble method could offer a rapid estimate of MMP, however, the MMP is determined by the visual operation, and thus the accuracy greatly depends on the operator [10,13]. Later the vanishing interfacial tension (VIT) technique was developed as an alternative experimental method to measure the MMP [14]. The VIT test is rapid and has a rigid definition of 'miscibility'. That is, the MMP is taken as a pressure at which the interfacial tension (IFT) between injected gas and crude oil becomes zero [14,15]. Orr et al [16] and Nobakht et al [17] have discussed the reliability and accuracy of the VIT test, showing good agreements with the slim tube test.
To predict the MMP by the VIT method, the IFT between injected gas and crude oil is indispensable to know. The pendant drop and the capillary rise methods are the most general experimental approaches to measure the gas/oil IFT. However, the pendant drop approach fails at low IFTs since it is difficult to form pendant drops [12]. The capillary rise approach needs to measure the density of two phases, the contact angle with capillary tube, and the accurate inner diameter, which requires high precision equipment and makes it complex, and this method has difficulty to measure the IFTs of CO 2 -live oil system since it is easy to form bubbles [7]. Besides experimental methods, some researchers [18][19][20][21] used parachor model to calculate IFTs, which requires the equation of state (e.g., Peng-Robinson equation of state [22] for vapor-liquid equilibria calculation). However, this method has a limitation that its parameters are adjusted by experimental data for a given reservoir.
Compared with experimental methods and the parachor model method, molecular modeling can be utilized to obtain IFTs of injected gas/oil interfaces and meanwhile reveal the microstructure and more detailed properties of injected gas/oil system, which help to understand the process of miscible-phase displacement. For instance, Boek and coworkers investigated the IFTs of CO 2 -brine system at various temperature, pressure and brine salinity by molecular dynamics (MD) simulation, and found that the IFT decreases with increasing pressure under any temperature condition but increases linearly with the molality of the salt solution [23]. de Lara et al explored the structural and dynamic properties of crude oil at the molecular level, and their calculations suggest an increase in fluid diffusivity of the oil phase with the IFT reduction [24]. More recently, the effects of methane on the water/octane interface were studied by molecular modeling, which showed that enlarging the ratio of methane could increase the interfacial roughness and interfacial thickness of the water/oil interface [25]. However, although researchers have studied the microstructure of the interface and the change of the IFT with MD simulation providing molecular information and insights, little work reported how to estimate the MMP of CO 2 injected into crude oil by MD simulation and understand the underlying mechanism.
In this study, we carried out MD simulations for various injected gas-crude oil systems, aiming to obtain the MMP in a rapid and economic way and provide molecular details and insights to help to understand the mechanism of some influence factors on the MMP. Specifically, to calculate the IFT of injected gas/oil interface, we used modified Lorentz-Berthelot rules to describe the interactions between injected gas and crude oil, which can greatly solve the mismatch from the original force fields and produce accurate IFTs [26]. Then, we investigated the effects of injected gas components, oil components and reservoir temperature on the MMP. We adopted the density distribution of each component, a new parameter (density ratio) to reflect the injected gas oleophilicity and the injected gas solubility in oil to explain their impacts. Finally, we conducted sensitivity analysis to evaluate the effects of those factors on the MMP.

Methods
As illustrated in figure 1, the MD simulation system consists of a slab of liquid oil phase in the middle and two slabs in gas phase on both sides of oil slab to form injected gas/oil interfaces. Injected gas considered in this study includes CO 2 , N 2 , methane, and their mixtures, and oil phase contains different types of alkanes and mixtures, depending on the system. To keep it simple, we abbreviate the alkanes involved in this work by their number of carbon atoms. For example, C 1 , C 6 and C 10 represent methane, n-hexane and n-decane, respectively. Specifically, for CO 2 +C 6 binary system, the oil slab contains 800 C 6 molecules while its two sides are occupied by different number of CO 2 molecules to produce different system pressure. For CO 2 +C 6 +C 10 ternary system, the total number of C 6 and C 10 in oil slab is 700 while their ratio is changed to explore the effects of the mole fraction of C 10 in oil slab on the MMP. For those systems with crude oil, the total number of all alkanes is 600, and detailed settings are listed in table S1. The force field for CO 2 is taken from [27], the force field for N 2 is taken from [28] and the NERD force field [29] is used for all alkanes since they exhibited more accuracy to describe the pure systems respectively.
All simulations were performed in the NVT ensemble using an MD package Gromacs, version 4.6.7 [30], with a box size of 5×5×35 nm 3 (figure 1). Periodic boundary conditions were used in three dimensions. The system temperature was maintained at a setting value with the Berendsen thermostat [31]. The electrostatic interactions were computed using the PME method [32]. An FFT grid spacing of 0.12 nm and cubic interpolation for charge distribution were taken to compute the electrostatic interactions in reciprocal space. A cutoff distance of 2.0 nm was used in the calculation of electrostatic interactions in the real space. The nonelectrostatic interactions were computed by direct summation with a cutoff length of 2.0 nm. The Leap-Frog algorithm was applied with 1 fs as time step to integrate Newton's equation of motion. The center of mass translation was removed every time step. Each simulation consists of an equilibration run of 10 ns and a production run of 10 ns. The IFT and pressure were calculated and stored every 0.1 ps using our developed onfly code. To ensure the simulation is well equilibrated within 10 ns, a test simulation containing 10 ns equilibration and 50 ns production was carried out. Results are displayed in figure S1, and show that the IFTs have notable fluctuation no matter how long the equilibration stage is, whereas the MMP fitted from the IFTs is of high accuracy. It tells us that plenty of IFT data points are needed for an MMP fitting (no less than eight IFT data points) and 10 ns is enough for equilibration to our system.
Since the Lorentz-Berthelot rules are the most common combination rule, it is adopted in the simulation of this work [33,34]. However, many researchers had reported that the standard Lorentz-Berthelot rules may misestimate the strength of mixture interactions, especially when the force fields were taken from different references [33,35,36]. The mixture interactions between injected gas and oil directly affected the IFTs. Therefore, we adopted a scaling parameter to revise the interactions as shown in equation (1).
This method exhibited a good accuracy to describe the gas/alkanes interfacial properties [26]. Here, a=0.9 was adopted for CO 2 -alkanes interactions, and similarly, a=0.85 was used for N 2 -alkanes interactions.
The IFT of injected gas/oil interface was calculated using the formulation of the Gibbs interfacial tension [24]. Two interfaces in simulation box were parallel to the XY plane (figure 1), thus the IFT was evaluated from the expression of pressure tensor as follows [37], where P N and P T represent the normal and tangential pressures, respectively, P αα (α=x, y, z) is the diagonal elements of the pressure tensor, and L z is the length of the simulation box in Z-direction.

Interfacial properties and MMP
According to the VIT method, the MMP is calculated by extrapolating the IFT at zero. Here, we conducted MD simulations for CO 2 +C 6 binary system, CO 2 +C 6 +C 10 ternary system and CO 2 +crude oil system, to obtain the MMP and provide molecular insights. The IFT between CO 2 and C 6 as a function of pressure at 80°C is shown in figure 2(a), which presents an approximately linear decrease with the increase of pressure in a wide range. Hence we take linear regression to get the MMP of this binary system and found it to be around 10.2 MPa, very close to the experimental measurement (10.8 MPa [38]), which demonstrates that our MD modeling has good feasibility and accuracy to estimate the MMP. A ternary system (CO 2 +C 6 +C 10 ) at 38°C was further studied, as figure 2(b) displays the MMP of this ternary system. When the mole fraction of C 10 in oil phase reaches 100%, the deviation of the simulation result from experimental data of Reamer et al [39] is about 0.6 MPa. Moreover, the MD-obtained MMP is found to increase from 6.5 to 8.6 MPa with the mole fraction of C 10 in oil phase rising. The effects of C 10 on the MMP is mainly due to its higher IFT of CO 2 /C 10 interface comparing with CO 2 /C 6 interface, thus the ratio of C 10 in oil phase increases the MMP of CO 2 +C 6 +C 10 ternary system.
Then Bakken oil with multiple components [7] was chosen as crude oil to validate our method for estimating the MMP. Oil components and their number of molecules in MD simulation are listed in table S1. Figure 3 exhibits the density profile of CO 2 -crude oil system. C 4 and longer alkanes have similar distribution along Z-direction, therefore, they are summed together as one component in figure 3. Specifically, figures 3(b) and (c) show the mass density for CO 2 -Bakken live oil system and CO 2 -Bakken dead oil system, respectively. It can be found that C 1 is denser in gas phase than in oil phase, which is similar to the distribution of CO 2 . Therefore, C 1 can be considered as an impure injection (CO 2 /C 1 mixture injection) rather than pure CO 2 injection. The C 1 in CO 2 will increase the MMP (see section 3.2.1), which accounts for a higher MMP for CO 2 -live oil system with presence of C 1 than CO 2 -dead oil system absent of C 1 .
MMPs of CO 2 -Bakken oil system predicted by MD simulation under different temperature are listed in table 1 (IFTs as a function of pressure are shown in figure S2). The MMP of CO 2 -Bakken live oil system determined by MD modeling agrees well with the slim tube test result at 129°C [7], and MMPs of CO 2 -Bakken dead oil system from this study are consistent with the VIT method data at two different temperatures [7]. Furthermore, it can be seen that the MMP increases with the temperature increasing for CO 2 -crude oil system, which has been reported by previous work [40][41][42]. MMPs of CO 2 -live oil system are higher than of CO 2 -dead oil system at the same temperature, coinciding with experimental results from Menouar et al [43], which may be due to the presence of C 1 in the former. Since as discussed above, the existence of C 1 will increase the CO 2 -crude oil MMP, the MMP is higher for CO 2 -live oil system than CO 2 -dead oil system. Compared with experimental data, MD-obtained MMPs have absolute relative deviations (ARD) within 5%, indicating that MD simulation could make acceptable estimation of the MMP for CO 2 -crude oil system. Therefore, considering its many advantages, such as efficiency, low cost, accessibility to a wide range of conditions and so on, MD simulation could be an alternative and competitive method for MMP estimation.

The influence factors of the MMP
After exhibiting the capability of MD simulation to estimate the MMP, we scrutinized how injected gas components, oil components and reservoir temperature could impact the MMP. CO 2 /C 1 mixture and CO 2 /N 2 mixture were taken as injected gas to reveal the gas component impacts. C 6 and C 19 were selected as representatives for intermediate components and heavy components in crude oil, respectively. The MMPs of N 2 -crude oil system and CO 2 -crude oil system at different temperature were calculated to observe the temperature influence. Figure 4 displays the effects of injected gas components on the CO 2 -crude oil MMP (as listed in table S1, the crude oil is a mixture of C 6 , C 10 , C 19 and C 30 ). The MMP of CO 2 /C 1 mixture with crude oil at various C 1 mole fractions in injected gas phase at 110°C (Figure 4(a)) shows that the MMP of this system increases nearly linearly with C 1 mole fraction rising, in agreement with the experimental results [7]. Figure 4(b) illustrates the MMP of CO 2 /N 2 mixture with crude oil as a function of N 2 mole fraction at 110°C, and demonstrates an increase in the MMP with increasing N 2 mole fraction. The MMP values of CO 2 /N 2 mixture with crude oil estimated by MD Figure 3. Density distributions along Z-direction. (a) Simulation snapshot of CO 2 -crude oil system, (b) density profiles of CO 2 -Bakken live oil system, (c) density profiles of CO 2 -Bakken dead oil system. Z/L is a normalized factor of the box in the Z-direction, and L=35 nm for both systems. The densities of C 1 and C 2 were enlarged by 10 fold for clarity. simulation are well in line with experimental data [7]. The MMP changes from 16.7 to 33.5 MPa as the mole fraction of N 2 rises from 0 to 0.5 while it just varies from 16.7 to 23.2 MPa for CO 2 /C 1 mixture, implying that N 2 has more significant impacts on the MMP of CO 2 -crude oil system than C1. It is worth noting that modified Lorentz-Berthelot rules play an important role in the MMP prediction. As shown in figure S3, the MMP and interfacial properties (such as mole fractions, IFTs, etc) deviate from experiment a lot, using standard Lorentz-Berthelot rules for N 2 -crude oil system, but they agree well with experiment when adopting modified Lorentz-Berthelot rules.

The effects of injected gas components on MMP
To further understand the effects of injected gas components on the MMP, we calculated the density ratio for a certain injected gas to quantitatively reflect its oleophilicity. The ratio was defined as the injected gas (CO 2 , N 2 , or C 1 ) density in gas phase divided by its density in oil phase. Obviously, the smaller the ratio was (i.e., at the same condition, there was more injected gas in oil phase), the more oleophilic this injected gas component was. From density profiles shown in figures 5(a), (b), it is found that such ratios for CO 2 , C 1 and N 2 are 1.45, 1.83 and 5.16, respectively, which suggests that CO 2 is more oleophilic than C 1 , and C 1 is more oleophilic than N 2 . This is why N 2 and C 1 will increase the MMP of CO 2 -crude oil system, and N 2 has stronger impacts on the MMP than C 1 . The different ratio for CO 2 , C 1 and N 2 could be attributed to their different interaction energy with crude oil. The gas-oil interaction energy (figure 5(c)) indicates that since CO 2 has strongest attractive interaction with crude oil, it is more oleophilic than N 2 and C 1 and thus has the smallest density ratio. Figure 6 exhibits the effects of crude oil components on MMPs. It can be found that C 6 and C 19 have opposite effects on CO 2 -crude oil MMP (the crude oil is a mixture of C 6 , C 10 , C 19 and C 30 as listed in table S2). The MMP decreases with the increase of C 6 mole fraction but increases with the C 19 concentration. Compared with C 19 , C 6 has a stronger impact on the MMP. The effects of C 6 and C 19 on the MMP are similar with the MMP influenced by C 6 and C 10 as shown in figure 2(b), that is, an oil component will change the MMP toward the MMP of the system of such a component and CO 2 when its mole fraction increases. Since the MMP of C 6 with CO 2 is 11.0 MPa which is lower than the one of crude oil (17.43 MPa) with CO 2 , C 6 reduces the MMP of crude oil with CO 2 , while C 19 increases the MMP.

The effects of oil components on MMP
The oleophilicity of CO 2 with an oil component also can explain the influence of oil components. As discussed above, we calculated the density ratio of CO 2 for different CO 2 +alkane (C 6 , C 10 , C 19 or C 30 ) binary systems ( figure S4). The ratio for CO 2 +C 6 binary system is 0.88, much smaller than that (1.45) for CO 2 with crude oil system; the ratio smaller than 1 represents that more CO 2 could be distributed in oil phase than in gas phase, implying strong solubility and oleophilicity of CO 2 in C 6 . Thus, C 6 is easier to be miscible with CO 2 , and can reduce the MMP of crude oil with CO 2 . Conversely, C 19 , showing a higher ratio for CO 2 +C 19 binary system than CO 2 with crude oil system, would increase the MMP of crude oil. The effects of a certain oil component on the MMP depend on the crude oil properties. Generally speaking, heavy oil components would increase while light oil components may decrease the CO 2 -crude oil MMP [44,45].

The effects of temperature on MMP
With exploring the temperature dependence, the MMP of CO 2 -crude oil system is revealed to increase with temperature, while the MMP of N 2 -crude oil system decreases during temperature increasing ( figure 7). This is caused by the different effects of temperature on the IFTs for these two systems. The IFT of CO 2 -crude oil system decreases with temperature at low pressure, while it rises with temperature at high pressures as shown in figure 7(a). In contrast to CO 2 -crude oil system, the IFT of N 2 -crude oil system reduces with increasing temperature at any pressure ( figure 7(b)). These phenomena have also been observed in experiments [40,41,46].
To further understand this different IFT trend between CO 2 and N 2 , the solubility in oil was evaluated by their density profiles, that is, gas density in oil phase is divided by oil density in oil phase. Previous work has reported that, for a gas/liquid interface, the IFT usually decreases with increased temperature or dissolved gas   [7]. The dashed line is just a guide to the eye. [47]. figure 8 shows the solubility of CO 2 and N 2 in oil as a function of temperature at various pressures, from which, it can be found that increasing temperature decreases the solubility of CO 2 but increases the solubility of N 2 in oil. Their different solubility can be directly observed and understood from their different density profiles as shown in figure 9. With temperature rising, the CO 2 density in oil phase decreases and meanwhile the oil density changes little, as a result, the CO 2 solubility in oil reduces with temperature. However, for N 2 -crude oil system, when temperature rising, the N 2 density in oil phase increases slightly while the oil density decreases due to the expansion of oil into N 2 phase, implying an increase in N 2 solubility with temperature. The abnormal solubility of N 2 in oil has also been reported in experiment [46,48]. Considering their different solubility trend, for CO 2 , the impacts of temperature and solubility on the IFT are opposite. At low pressure, the impact of temperature on the IFT is dominant since the solubility is small and changes little, that is, with the near same amount of CO 2 and oil at gas-liquid interfaces, high thermal motion would facilitate the CO 2 -oil mixing and then to reduce IFT. However, at high pressure, the reduction of CO 2 solubility becomes dominant, indicating an increase of IFTs that would impede the CO 2 -oil mixing. For N 2 , both temperature and dissolved gas have the similar effects on the IFT, that is, decrease the IFT of N 2 -crude oil system.

Sensitivity analysis
Sensitivity analysis is a widely used approach to evaluate the effects of related parameters (in this work, those would include reservoir temperature, C 6 and C 19 mole fraction in oil phase, C 1 and N 2 mole fraction in injected gas phase) on the MMP of CO 2 -crude oil system [49]. The correlation coefficient is usually calculated for sensitivity analysis, which can be evaluated using the following formula [6, 50], where p k i , and p k are the ith value and the average value of the kth influence parameter, respectively; MMP i and MMP denote the corresponding MMP and average value, respectively.
The correlation coefficient can assess the influence degree of each parameter on the MMP. Higher absolute value of the correlation coefficient suggests greater impact on the MMP, and positive value means this parameter facilitates to increase the MMP.
Results are illustrated in figure 10, which reveal that all parameters except C 6 mole fraction will increase the CO 2 -crude oil MMP. N 2 mole fraction and temperature have considerable effects on the MMP since they have higher absolute value. The impurities in the injected CO 2 gas (i.e., N 2 or C 1 ) bring a notable increase in the MMP, and N 2 mole fraction has more significant influence on the MMP than C 1 mole fraction (see figure 4). This can attribute to the fact that N 2 and C 1 have lower oleophilicity than CO 2 as discussed in section 3.2.1. As to oil components, C 6 mole fraction will reduce the MMP while C 19 increases it, and C 6 has stronger effects on the  . Sensitivity analysis to find the impacts of those parameters on the CO 2 -crude oil MMP. The parameters include injected gas components (C 1 and N 2 ), oil components (C 6 and C 19 ), and temperature (T).
MMP than C 19 which also can be found in the previous analysis (see section 3.2.2). This implies that heavy oil components trend to increase the MMP while intermediate oil components decrease it.

Conclusions
We conducted MD simulations to estimate the MMP of CO 2 -crude oil system based on calculations and analyses of IFTs between injected gas and crude oil. In order to obtain accurate IFTs, modified Lorentz-Berthelot rules were developed to properly describe the interactions between injected gas and oil. The MMPs of CO 2 +alkane binary system and more complex system like CO 2 -crude oil system predicted by MD modeling are in good agreement with experimental results, which demonstrates that MD simulation has good feasibility and accuracy to estimate the MMP, making it an alternative, and competitive method for the MMP estimation.
The effects of injected gas components, oil components and reservoir temperature on the MMP were investigated and results revealed that: (1) the ratios of C 6 and C 19 have opposite effects on CO 2 -crude oil MMP, that is, the MMP increases with the decrease of C 6 mole fraction or the increase of C 19 mole fraction; (2) C 1 and N 2 will increase the MMP of CO 2 -crude oil system, and N 2 has stronger impact on the MMP than C 1 ; (3) The MMP of CO 2 -crude oil system increases with temperature, while the MMP of N 2 -crude oil system decreases with temperature increasing. This is due to their different solubility in oil at various temperatures, leading to different IFT changes with temperature. Finally, by the sensitivity analysis, we can find that the concentration of N 2 in CO 2 and reservoir temperature have considerable effects on the CO 2 MMP. Mole fraction of C 6 will decrease the MMP while the others increase the MMP.
This work potentially renders a simple and rapid way (i.e., MD simulation, and a cluster containing 20×24 CPU cores can obtain about one MMP value per 8 h) compared with days-long experimental methods [11] to estimate the MMP, and the discussion of the impacts of some parameters (T, gas components, oil components) on the MMP can help to understand the process of miscible-phase displacement. The VIT method does not consider the effect of rocks, however, the effect of nanopores sometimes is not negligible, especially for unconventional reservoirs. Thus, the MMP in nanopores would be important and need further study.