Bulk and interfacial properties of decane in the presence of carbon dioxide, methane, and their mixture

Molecular dynamics simulations were performed to study the bulk and interfacial properties of methane + n-decane, carbon dioxide + n-decane, and methane + carbon dioxide + n-decane systems under geological conditions. In addition, theoretical calculations using the predictive Peng-Robinson equation of state and density gradient theory are carried out to compare with the simulation data. A key finding is the preferential dissolution in the decane-rich phase and adsorption at the interface for carbon dioxide from the methane/carbon dioxide mixture. In general, both the gas solubility and the swelling factor increase with increasing pressure and decreasing temperature. Interestingly, the methane solubility and the swelling of the methane + n-decane system are not strongly influenced by temperature. Our results also show that the presence of methane increases the interfacial tension (IFT) of the carbon dioxide + n-decane system. Typically, the IFT of the studied systems decreases with increasing pressure and temperature. The relatively higher surface excess of the carbon dioxide + n-decane system results in a steeper decrease in its IFT as a function of pressure. Such systematic investigations may help to understand the behavior of the carbon dioxide-oil system in the presence of impurities such as methane for the design and operation of carbon capture and storage and enhanced oil recovery processes.

oil 24 . Traditionally, natural sources of CO 2 are used for about 90% of EOR production in the U.S 23 . However, EOR using anthropogenic CO 2 emissions would be necessary to achieve the desired climate benefits. CO 2 captured from power plants and industrial sources may contain impurities [27][28][29] . The cost of purification of gas mixtures considerably increases for obtaining a product of high degree of purity. Therefore, it is important to the study the behavior of the CO 2 -oil system in the presence of impurities such as methane for the design and operation of CCS and EOR processes.
Molecular simulations provide valuable insights into the bulk and interfacial properties of various systems [14][15][16][17][58][59][60][61][62][63] . In this work, we perform molecular dynamics (MD) simulations to understand the properties of CH 4 + n-decane, CO 2 + n-decane, and CH 4 + CO 2 + n-decane systems over a broad range of temperature (313-442 K) and pressure (up to about 300 bar) relevant to geological processes. Given the general lack of experimental results for these systems, especially the ternary case, theoretical modeling is employed to complement the simulation data.

Methods
Simulation details. The MD simulation was carried out using the GROMACS code 64 . The simulation model and method used here are similar to the ones reported in our previous work 14,63 . The transferable potentials for phase equilibria force field 65 was used to model n-decane, methane, and CO 2 . The simulations were initiated by placing 400 decane molecules and up to 2048 CH 4 /CO 2 molecules in each simulation box. In directions parallel to the interface, each cell was about 50 × 50 Å, and periodic boundary conditions were used in all directions (Fig. 1). These values were chosen to ensure reasonable bulk phases with enough molecules and also to reduce the truncation and system size effects in the calculations of the bulk and interfacial properties 39,46,[52][53][54][55][56] . These systems are energy-minimized using the steepest descent algorithm and subjected to equilibration MD runs of 10 ns in the NPT ensemble. During the NPT runs, the volume changes are obtained by adjusting only the cell size in the z direction (L z ) normal to the interface. After that we carried out production runs of 6 ns duration in the NVE ensemble. The leap-frog algorithm 66 with a time step of 1 fs was used to integrate the equations of motion. The linear constraint solver algorithm 67 was implemented to constrain the bond lengths and angles. Virtual sites were used to maintain the rigid structure of CO 2 68 . The particle mesh Ewald method was used to calculate the long-range electrostatic interactions. The Nosé-Hoover thermostat and Parrinello-Rahman barostat, with a relaxation time of 2.0 ps for both, were used to regulate the temperature and pressure, respectively. Each MD run was repeated three times and the average value was taken.
The IFT was calculated as follows 52,63 : z zz x x y y where P xx , P yy , and P zz are the three diagonal components of the pressure tensor. Typically, for the estimation of various atomic density profiles, a bin width of 1 Å is used.
theoretical details. The volume translated predictive Peng-Robinson 1978 (VT-PPR78) 40,44,69 equation of state (EoS) with the usual van der Waals one-fluid mixing rules and a linear mixing rule for the volume correction is used to estimate bulk properties. The VT-PPR78 model calculates the binary interaction parameter k ij for each temperature using a group contribution method. Additionally, the interfacial properties were calculated by coupling the VT-PPR78 EoS with the density gradient theory (DGT). In brief, the Helmholtz free energy of the planar interface is given by 70,71 Figure 1. Equilibrium snapshot of the CH 4 + CO 2 + n-decane system at 343 K and 60 bar.
where A is the area of the planar interface, f 0 represents the Helmholtz free energy density of the homogeneous fluid at the local density n, dn i /dz is the local gradient in density of a given component i, and c ij is the cross influence parameter defined as follows: where c ii and c jj denote the pure component influence parameters, and β ij represents the binary interaction coefficient. In this study, the pure component influence parameter is computed based on the correlation developed by Miqueu et al. 44 . The minimization of the Helmholtz free energy leads to the following Euler-Lagrange equations: for , 1, , , , μ i denotes the chemical potential of component i in the equilibrium bulk phase and N c represents the total number of components. To obtain the density profiles across the interface, the above set of equations must be solved together with the following boundary conditions: where n i I and n i II are the equilibrium densities of component i in the coexisting bulk phases. As stated above, the bulk properties are calculted using the VT-PPR78 EoS. Using the geometric mixing rule (β ij =0) reduces the system of N c differential equations (Eq. (4)) to the system of (N c − 1) algebraic equations: ) ( ( ( )) ); 1, , , and ref These equations could compute the densities of all except one component and the density of the latter is monotonic function of z over the whole interface 40,44 and used as a reference variable. Once the density profiles are known, the IFT can be calculated using the following equation: at the local density and its value in the bulk phase Ω B = −P. Furthermore, the z-coordinate for the density profile can be computed using the following equation:

Results
Atomic density profiles. Both MD simulations and DFT calculations were used to predict the atomic density profiles that are not always easily accessible to experiments. For instance, Fig. 2 displays the density profiles of various species as obtained from the MD simulations (solid lines) and the corresponding theoretical data (dashed lines) for CH 4 + n-decane, CO 2 + n-decane, and CH 4 + CO 2 + n-decane (equimolar CH 4 /CO 2 mixture) systems at 343 K. We see that the simulation results are in good qualitative agreement with the theoretical data.
Note that among the different species, the maximum difference (differ by a factor of up to ~2) between the simulated and calculated densities is obtained for CO 2 , especially at high pressures. In all cases, the density profiles show local enrichment of gas molecules at the interface, but no such behavior is seen for n-decane in agreement with the previous studies [38][39][40][41][43][44][45][46][47][48][49][50][51][52][53] . That is, a positive surface activity is seen for both methane and carbon dioxide (dn i /dz = 0; d 2 n i /dz 2 < 0 in the interfacial region). Whereas the density profile of decane across the interface varies monotonically and does not show surface activity. Also, the higher methane and/or carbon dioxide contents result in a broadening of the interface. As seen from the adsorption peaks in Fig. 2, the local accumulations of methane and carbon dioxide, in general, increase with pressure. However, at high pressures, the local enrichment of gas molecules decreases with pressure. It is known that the local enrichment of the interface in methane and carbon dioxide typically decreases with temperature 38-40,43-49,52,53 . The density of CO 2 is usually higher in the liquid phase (decane-rich) when compared to that in the CO 2 -rich phase, whereas an opposite trend is found for methane. Our results also show the preferential dissolution in the decane-rich phase and adsorption at the interface for carbon dioxide from the CH 4 /CO 2 mixture.
Bulk properties. Figure 3 shows the bulk densities as obtained from the MD simulations (open symbols) and the corresponding theoretical data (solid lines) for CH 4 + n-decane, CO 2 + n-decane, and CH 4 + CO 2 + n-decane (equimolar CH 4 /CO 2 mixture) systems under geological conditions. These results are also compared with the corresponding experimental data 40 (closed symbols). Other experimental data 30,32,33,35,39 are similar to the ones plotted here and not shown for clarity. Note that in the absence of any experimental data, the simulation results of the CH 4 + CO 2 + n-decane system are only compared with the theoretical predictions. We see that both our theoretical and simulation results are in good agreement with the experimental data. For example, the overall absolute average deviation between theoretical and experimental density values is less than about 2% in the decane-rich phase, while it is less than about 3% in the CH 4 -rich or CO 2 -rich phase. In all cases, these total densities generally decrease with temperature in both phases. Also, these densities generally decrease with pressure in the decane-rich phase and increase with pressure in the CH 4 -rich and/or CO 2 -rich phases. However, for example, in the CO 2 + n-decane system, the density of the decane-rich phase increases with pressure, especially at low temperatures and pressures. These different scenarios may take place depending on the gas solubility and swelling in the decane-rich phase. The knowledge of the species mole fractions and swelling effects plays an important role in the use of carbon dioxide for the recovery of hydrocarbon fluids from underground reservoirs 30,32,40,53,56 . Figure 4 shows the bulk mole fractions of methane and carbon dioxide in the decane-rich phase as obtained from the MD simulations (open symbols) and the corresponding theoretical data (solid lines) for all studied systems. These results are also compared with the experimental data 30,32 (closed symbols). Both our theoretical and simulation results are again consistent with the experimental data. In general, various mole fractions for all studied systems increase with increasing pressure and decreasing temperature. Notably, the mole fraction of methane in these systems is not significantly affected by temperature. Overall, the amount of n-decane in the CH 4 -rich and/or CO 2 -rich phases was estimated to be relatively small (see, e.g., Fig. 2).
Furthermore, Fig. 5 shows the swelling factor for the decane-rich phase as obtained from the MD simulations (open symbols) and the corresponding theoretical data (solid lines) for all studied systems. The swelling factor is defined as ratio of the volume of the saturated decane-rich system to the volume of the decane alone 53,56 . Though not shown, both theoretical and simulation results for densities of the pure n-decane system at the studied conditions are in good quantitative agreement with the NIST data 72 . In general, the swelling factor for all the studied systems increases with increasing pressure and decreasing temperature, which is consistent with the above solubility data. Also, similar to the solubility vs temperature behavior, the swelling of the CH 4 + n-decane system is not significantly affected by temperature. It is known that the interaction between CO 2 and alkane molecules plays a key role in the solubility and swelling properties of the alkane-CO 2 system 53,56 . For example, Fig. 6 shows the radial distribution functions between decane (carbon site) and CH 4 /CO 2 in the bulk liquid phase of the CH 4 + CO 2 + n-decane system at 343 K and 60 bar. These results were obtained using an NPT MD simulation of the corresponding bulk liquid phase alone. From this figure, it can be seen that CO 2 molecules are positioned relatively close to decane. Furthermore, for this system, we find that the interaction energies of decane-CO 2 and decane-CH 4 are about −2428.5 and −930.9 kJ/mol, respectively. All these results indicate the preferential interaction of decane with CO 2 relative to CH 4 .
The conformation of the n-decane molecule may be characterized by its radius of gyration R g 54,55,61 . The R g is defined as the root-mean-square distance of all atoms from their common center of mass. The R g was calculated www.nature.com/scientificreports www.nature.com/scientificreports/ from the MD simulations in order to check the compactness of the conformations. Our calculations find that the R g values decrease with increasing temperature and are in the range of 0.32-0.34 nm for all systems. This is because as the temperature increases, the effective monomer-monomer attraction becomes increasingly important, causing the chains to shrink 54 . The values of R g are in good agreement with the previous reported results 54,55 . We also find that the R g is much less sensitive to changes in pressure and to the presence of methane and CO 2 .
interfacial properties. Figure 7 shows the IFTs as obtained from the MD simulations (open symbols) and the corresponding theoretical data (solid lines) for CH 4 + n-decane, CO 2 + n-decane, and CH 4 + CO 2 + n-decane (equimolar CH 4 /CO 2 mixture) systems under geological conditions. These results are also compared with the corresponding experimental data 40 (closed symbols). Other experimental data 31,[33][34][35]37,39,41 are similar to the ones plotted here and not shown for clarity. Note that in the absence of any experimental data, the simulation results of the CH 4 + CO 2 + n-decane system are only compared with the theoretical predictions. In the temperature range (313-442 K) considered here, the surface tension of pure n-decane decreases from about 22 to 11 mN/m 72 . It can be seen that both our theoretical and simulation results are in good agreement with the experimental estimates. For example, the overall absolute average deviation between theoretical and experimental IFT values is less than about 9%. These results show that, at a given temperature and pressure, the presence of CH 4 increases the IFT of the CO 2 + n-decane system. The IFT of CH 4 + n-decane, CO 2 + n-decane, and CH 4 + CO 2 + n-decane systems decreases with pressure. In general, a similar trend is observed for the temperature dependence of the IFT. However at high pressures, for example, the IFT of the CO 2 + n-decane system increases with increasing temperature. This can be described 40 by the fact that the local enrichment of gas molecules at the interface (see, e.g., Fig. 2) is more pronounced for carbon dioxide when compared to that of methane, especially at low temperatures. The correlation between the IFT and the local enrichment of the interface in methane and carbon dioxide will be discussed in detail below. Notably, theoretical calculations based on other equations of state 38-41,43-51 , such as the statistical associating fluid theory, can also give accurate estimates of the bulk and interfacial properties of these www.nature.com/scientificreports www.nature.com/scientificreports/ binary systems. One of the main properties governing the design and operation of CCS and EOR processes is the miscibility of the injected CO 2 with oil 24,36 . The minimum miscibility pressure is taken as a pressure at which the IFT becomes zero 36 . Due to the fact that the IFT of the CO 2 + n-decane system decreases steeply with increasing pressure, pure CO 2 is completely miscible with decane from relatively low pressures. In the presence of methane, the CO 2 + decane system seems to indicate a smaller miscibility region with respect to pressure.

Discussion
The results obtained in this study allow one to quantify the the local enrichment of the interface in methane and carbon dioxide. In this context, it is important to notice that the IFT is related to the surface excess Γ i as given by the Gibbs adsorption equation 40,63 :  4 and CO 2 in the decane-rich phase: (a) CH 4 in CH 4 + n-decane, (b) CO 2 in CO 2 + n-decane, (c) CH 4 in CH 4 + CO 2 + n-decane, and (d) CO 2 in CH 4 + CO 2 + n-decane systems. The symbols are the same as in Fig. 3. Note that here the experimental data were taken from Reamer et al. 30,32 .
The surface excess of a component i with respect to decane is given by Here, I represents the CH 4 -rich and/or CO 2 -rich phases, and II the decane-rich phase. Figure 8 shows the surface excesses of CH 4 and CO 2 as obtained from the MD simulations (open symbols) and the corresponding theoretical data (solid lines) for all studied systems. These results are also compared with the experimental Figure 6. Radial distribution functions between decane and CH 4 /CO 2 in the bulk liquid phase of the CH 4 + CO 2 + n-decane system at 343 K and 60 bar. www.nature.com/scientificreports www.nature.com/scientificreports/ data 35,40 (closed symbols). In all cases, we see a nonmonotonic variation of the surface excess with pressure. The relatively higher surface excess of the CO 2 + n-decane system results in a steeper decrease in its IFT as a function of pressure (see also Eq. (10) and Fig. 7). Overall, the surface excesses decrease with increasing temperature and the results of the ternary system fall between those of the binary systems.
To conclude, we have studied using MD simulations the bulk and interfacial properties of CH 4 + n-decane, CO 2 + n-decane, and CH 4 + CO 2 + n-decane systems under geological conditions. Theoretical modeling using the VT-PPR78 EoS combined with DGT was employed to complement the simulation data. The simulation results of the atomic density profiles and the IFT values are in good agreement with both the theoretical and experimental data. The results of our study show that, in general, the local accumulations of methane and carbon dioxide increase with increasing pressure and decreasing temperature. At high pressures, however, the local enrichment of gas molecules decreases with pressure. An important result is the preferential dissolution in the decane-rich phase and adsorption at the interface for carbon dioxide from the CH 4 /CO 2 mixture. The gas solubility and swelling in the decane-rich phase can play a key role in determining the behavior of the bulk densities as a function of temperature and pressure. Typically, both the gas solubility and the swelling factor for all the studied systems increase with increasing pressure and decreasing temperature. Notably, the mole fraction of methane and the swelling of the CH 4 + n-decane system are not significantly affected by temperature. Furthermore, we find that the presence of CH 4 increases the IFT of the CO 2 + n-decane system. The IFT of the investigated systems decreases with pressure. In general, a similar trend was observed for the temperature dependence of the IFT. However, at high pressures the IFT of, e.g., the CO 2 + n-decane system increases with temperature. There is a nonmonotonic variation of the surface excess with pressure and it decreases with temperature. The surface excess is more pronounced for carbon dioxide when compared to that of methane in the studied systems. Because of this fact, the IFT decreases steeply with increasing pressure for systems containing CO 2 .