How well can Damped Shifted Force Monte Carlo predict vapor-liquid equilibria for natural gas systems?

The oil and gas industry has faced new challenges, especially since discovering and exploring Brazilian pre-salt reserves. The current dehydration processes, imperative for gas processing and distribution, result in relevant solvent losses or are too expensive, space-demanding, or both. Molecular simulations have been among the preferred approaches in the last few years. Although computationally demanding, their versatility, predictive character, and accuracy make them an asset for many tasks. As the required time can yet be a hindrance when working with molecular simulations, the Damped Shifted Force (DSF) method can be an exciting option for a faster calculation of the electrostatic interactions, in contrast with the more popular Ewald methods, more rigorous and time-consuming. The present work presents Monte Carlo simulations using the DSF method to predict the vapor-liquid equilibria of binary mixtures containing some of the natural gas main components: CH 4 , CO 2 , and H 2 O. The results were compared to experimental data and other Monte Carlo simulations found in the literature to check whether the DSF method is a reliable alternative for natural gas systems. The comparison showed that DSF simulations had comparable accuracy to other Monte Carlo simulations (including Ewald ones) and provided predictions close to experimental data.


Introduction
Natural gas has been a signi cant fuel for decades, with an essential role in the energetic and economic scenario, with consistent increases in global consumption (BP 2021).The current concern towards developing more environmental-friendly processes puts natural gas in an even more highlighted position because it is cleaner than other fossil fuels, like petroleum and coal (Santos et al. 2021).
Removing contaminants, like water, carbon dioxide, and hydrogen sul de, is mandatory for natural gas distribution, both for reaching commercial speci cations and because these impurities favor corrosion, hydrate formation, and heat capacity loss (Santos et al. 2017;Zou et al. 2013;Gonfa et al. 2015).
The discovery of Brazilian pre-salt reserves promoted the rise of new challenges in the oil and gas industry.Particularly for natural gas production, the high gas-oil ratio and high carbon dioxide content are noteworthy (Aimoli et al. 2020).The traditional dehydration methods need to be revised in this scenario, given that they present relevant drawbacks, especially for offshore applications.For instance, glycol absorption results in a relevant solvent loss and foaming (Heym et al. 2010;Yu et al. 2021).In contrast, adsorption with zeolites struggles with the need for big spaces, and dehydration by condensation might require external cooling and hydrate inhibitors, leading to additional regeneration structures (Netusil and Ditl 2011).In the last decade, supersonic separation received extra attention.Although it provides an e cient, compact separation without using solvents or reagents, it promotes a substantial gas pressure drop (Yang et al. 2014; Shooshtari and Shahsavand 2017; Niknam et al. 2018).Finally, membrane separation is another alternative, proposed mainly for submarine gas treatment due to its absence of mobile parts, compact character, and low frequency of maintenance.However, the possibility of losing up to 1% of the methane is an important drawback compared to the other methods, especially for mild operation conditions (He et al. 2020;Ahmadi et al. 2021).
In this context, the proposition of novel solutions is desirable and has been pursued in this eld.For instance, glycol replacement in the absorption process with ionic liquids has been reported in the past years (Gonfa et al. 2015;Yu et al. 2021;Yu et al. 2017).The myriad of possible solvents and their mixtures and the necessity of knowledge about their properties make developing and applying predictive methods imperative.Among the predictive methods, molecular simulations draw some attention due to their versatility, the transferability of its potentials, and overall good performance.The success and consequent popularity of the molecular simulation methods resulted in many options of software and packages with diverse approaches, from more speci c to more general applications.Both commercial and free implementations are available.Shah et al. [16] provided an interesting summary of the most popular molecular simulation tools.
Among molecular simulation methods, Monte Carlo simulations are often used for studying phase equilibria.The major drawback of molecular simulation approaches is related to their high computational demand.The method used for calculating the electrostatic interactions between atoms in the simulation box is a critical factor that signi cantly in uences the simulation times.The preferred methods are the Ewald summation and their modi cations.These methods, although rigorous, are associated to computational expenses that can be a hindrance if the system is complex or if the computational resources are restricted.Fennell e Gezelter (2006) showed that the Shifted Force and Shifted Potential methods can be an alternative to Ewald methods, requiring less simulation time and comparable performance.Kolafa et al. (2008) pointed at the purpose of the prediction as a signi cant factor in deciding the electrostatic interactions calculation method.The authors point to the force eld as the most signi cant source of error.Hence simulations trying to reproduce experimental data should be able to use Shifted Force reliably and Shifted Potential methods.In opposition, simulations that develop new theories or force elds should rely on more accurate electrostatic interaction calculation methods.Among these methods, the Damped Shifted Force method (DSF) has shown satisfactory results, considering both performance and simulation time, even for reasonably complex systems In the present work, the vapor-liquid equilibria of binary systems of crucial natural gas constituents (methane, water, and carbon dioxide) were simulated through Monte Carlo simulations with the DSF as the electrostatic interactions calculation method.These results were compared to other approaches and experimental data to evaluate how good DSF simulations can predict phase equilibria for methane and its primary contaminants in a natural gas production process.

Simulation details
The systems were simulated in the CASSANDRA package (Shah et al. 2017) using the Gibbs Ensemble.
For all systems, the damping factor was set to 0.2 A − 1 , given that several reports points to that value as an optimal or close to optimal value (McCann and Acevedo 2013; Alcantara et al. 2020;Shi and Maginn 2008;Mei et al. 2016;Zahn et al. 2002).The cutoff radius was set to 12 A, as it can be a good compromise between time and performance even for more challenging systems (McCann and Acevedo 2013), for both van der Waals and electrostatic interactions, with analytic tail correction for the former.
Equilibration runs were done before the production runs.For an equilibration to be considered successful, it was decided it had to meet two criteria: each equilibration had to go through at least 100 million Monte Carlo steps and no systematic oscillations were detectable in the long term.Obviously, due to stochastic nature of the method, random oscillations will always be visible.For systematic oscillations, the equilibration was continued until this behavior was not detectable anymore.An example of equilibration with no long-term property variation can be seen in Fig. 1.After the equilibration was nished, production runs composed by 100 million Monte Carlo steps were done using the nal con guration of the equilibration run as the starting con guration of the production.The macroscopic properties of the system were obtained averaging the properties of the micro-states in the production.
The probability of each Monte Carlo movement was set as follows: 45% for translation, 45% for rotation, 9.5% for molecule swap between boxes and 0.5% for volume change.Such values were de ned after preliminary tests to achieve a compromise between simulation time and prediction quality.Three independent simulations were conducted for each condition to evaluate the related uncertainties.
The simulations generated in this work were compared with other simulations, experimental data and equations of state that are commonly available in commercial simulators (AspenTech 2017).The equations of state selection was done through a series of preliminary tests comparing their predictions with experimental data and is not the focus of this work.

Binary system CH 4 + CO 2
The selected force eld for carbon dioxide was the TraPPE (Potoff and Siepmann 2001), whereas the TraPPE-UA (Martin and Siepmann 1998) was used for methane.Aimoli et al. ( 2014) studied the prediction with those and other potentials and showed both TraPPE and TraPPE-UA can provide reliable predictions for many properties for pure carbon dioxide and methane, respectively, even for pressures up to 100 MPa.
For this system, the equilibria were simulated using the NVT Gibbs Ensemble.Although the NPT Gibbs Ensemble is more common for predictions for systems containing mixtures, this type of prediction can be done with both ensembles (Panagiotopoulos 1992;Ramdin et al. 2016).A total of 1000 molecules was used.

Binary system CO 2 + H 2 O
For water molecules, the force eld was the SPC/E (Berendsen et al. 1987), whereas CO 2 was still represented by TraPPE.The simulation was done in the NPT Gibbs Ensemble, with a total of 1000 molecules (700 water molecules and 300 CO 2 , similar to what was proposed by Vorholz et al. (2000)).

Binary system CH 4 + H 2 O
The TraPPE-UA and SPC/E force elds were used to represent methane and water, respectively, in the NPT Gibbs Ensemble.However, simulations with 1000 total molecules showed great di culty in providing good predictions, due to the very low solubility of each substance into the other, for both phases, as it can be seen in the data reported by Frost et al. (2014).Thus, 2000 molecules were used for this system.The increase in the total number of molecules brought the need to also change the equilibration procedure, given that larger systems take more Monte Carlo steps to equilibrate.For this system, the minimal amount of equilibration steps was increased to 200 million.

Binary system CH 4 + CO 2
The CH 4 + CO 2 vapor-liquid equilibrium curve was rst simulated at 250 K and the results were compared to experimental data provided by Wei et al. (1995).Additionally, the values predicted by Do et al. (2010) with Monte Carlo simulations were also used for comparison.The results are presented in Fig. 2.
It is important to say Do et al. chose an explicit hydrogen force eld to represent methane, the TraPPE-EH, and used the EPM force eld for CO 2 .They also used a spherical cutoff of 12 A with tail correction even for the electrostatic contribution.This approach towards the electrostatic interactions reduces the computational time, allowing the use of a more rigorous force eld, like the TraPPE-EH.However, this simple truncation can introduce relevant prediction errors associated with incorrect representation of the electrostatic interactions.It is convenient to remember that the electrostatic interactions decay slower with the distance than Lennard-Jones interactions, what makes a simple truncation a riskier approach for electrostatic interactions.
It is possible to see that the DSF predictions were superior to those reported by Do et al. for the dew point section of the equilibrium curve.One could expect the rigorous force eld makes up for the less rigorous electrostatic calculation method, but in this region of the curve, it is not what is observed.The DSF method, more rigorous than a simple truncation, performed better even with a united-atom force eld.Both approaches overestimated slightly the methane content for that region.
For the bubble point curve, the predictions from Do et al. were superior, being very close to the experimental data.DSF simulations, although with inferior performance, could correctly describe the equilibrium.Considering the whole curve, the GERG-2008 equation of state (Kunz and Wagner 2012) provided the best description of the equilibrium.It is important to remember that this equation was parameterized considering 20 main components present in natural gas streams, therefore its predictivity is restricted when introducing other components, limiting the search for new solutions in the specialized industry.
The CH 4 + CO 2 binary equilibrium was also simulated at 270 K and the results were compared to the experimental data reported by Wei et al. (1995) and to Monte Carlo simulations provided by Ramdin et al. (2016).They are shown in Fig. 3.In this case, the authors used the same force elds as the present work, but the electrostatic calculations were done by the Ewald method.
For this temperature, the comparison is more direct, given that the same force elds were used.The simulations of the present work and those by Ramdin et al. differ mostly by the electrostatic interactions calculation method and by the cutoff radius, as Ramdin et al. used the Ewald method and 14 A, respectively.For both the bubble point and dew point curves, DSF simulations provided better predictions.This result may sound unexpected, give the more rigorous character of the Ewald method.Some factors may explain this, like the cutoff radius.Preliminary simulations for this work using the same methodology showed better results for a cutoff radius of 12 A rather than 14 A, which may point at a systematic overestimation of the methane content with the increase of the cutoff radius.
On the other hand, Monte Carlo methods are inherently in uenced by some degree of stochasticity.It can be seen that DSF predictions were not vastly superior, thus when the uncertainties are considered, it is reasonable to say the predictions are approximately equivalent.Unfortunately, Ramdin et al. don't provide such uncertainties.The uncertainties of the present work, as well as the values provided by the simulations, are provided in Supplementary Information.
Finally, it is possible that small differences are found when using different simulators, as shown by Gowers et al. (2017).Ramdin et al. used the RASPA simulator (Dubbeldam et al. 2017).Anyway, the DSF simulations were able to predict competently the vapor-liquid equilibria of the CH4 + CO2 binary system, being quantitatively comparable to experimental data, equation of state and Monte Carlo simulations with Ewald method.
Another property that is important for the design of new processes is the density.The quality of the predictions for the density was also investigated and is represented in Figs. 4 and  Unfortunately, they didn't provide the exact values found in experiments, thus a more in-depth, accurate comparison is di cult.Like for the description of the equilibria compositions, the DSF simulations could describe the property correctly, with quality comparable to more rigorous methods, and the data, as well as the uncertainties, can be found in Supplementary Information.

Binary system CO 2 + H 2 O
In the case of the CO 2 + H 2 O binary, the vapor-liquid equilibrium was simulated at 348.15K and compared to experimental data (Zawisza and Boguslawa 1981;Wiebe 1941;Coan and King Jr. 1971; Wiebe and Gaddy 1939), as well as Monte Carlo simulations by Vorholz et al. (2000) and the predictive Soave-Redlich-Kwong and Schwartzentruber-Renon equations of state (AspenTech 2017).The results are presented in Figs. 6 and 7. Vorholz simulated this system using the EPM and SPC potentials for carbon dioxide and water, respectively.
Like what was seen for the rst system at 250 K, the DSF simulations, despite not giving the most accurate predictions among all analyzed, were capable of both correctly representing the behavior of the system and providing predictions quantitatively comparable to Ewald simulations and equation of state.
There is, though, a slight underestimation of the CO 2 mole fraction that increases with increasing pressure, what can also be observed for the simulation data by Vorholz et al (2000).
For the study of these data, it is important to point certain controversy regarding the experimental values.Aasen et al. (2017) suggested that those bubble point molar fractions reported by Wiebe and Gaddy (1939) are, actually, slightly under the real values.They point to this trend through the evaluation and comparison of other literature reports.On the other hand, Diamond and Akin ev (2003), doing an analysis of data quality from several references, attributed the highest level of con dence to the values published by Wiebe and Gaddy.Another relevant point is that the study by Aasen et al. comprised only the temperatures of 298.1 and 323 K and the alleged deviations seem to decrease with increasing temperature.Hence, it is expected that the deviations, already small even for lower temperatures, are even smaller for 348.15K, if they really exist.
For the CO 2 -rich region, the experimental data present a notable behavior of stabilization close to 0.99 CO 2 mole fraction with increasing pressure.For even higher pressures, even a decrease in this fraction can be observed, as shown by Aasen et al.Besides providing predictions with deviations comparable to other approaches, the DSF simulations were capable of correctly representing the mole fraction stabilization behavior, what doesn't happen when using the equation of state, for instance.
The density of each phase was also simulated and is shown as a function of pressure in Fig. 8.The predictions by other approaches are provided as well.What draws attention is the severe disagreement between the PSRK equation of state and the simulations for the liquid phase density.Though, the incapability of reliably predicting the liquid phase density is a well-known aw of the SRK equation of state and its modi cations (Ashour et al. 2011).Because of that, the Peng-Robinson (PR) equation of state was also included in the comparison.The exact values of the simulations and the uncertainties can be found in the Supplementary Information.

Binary system CH 4 + H 2 O
The prediction of CH 4 + H 2 O equilibrium turned out to be especially tricky, because the non-ideality of the system results in extremely low mutual solubilities.This fact resulted in the need of increasing the number of the molecules in the simulation boxes to avoid nite size effects, what, ultimately, led to simulations multiple times more computationally demanding.The results are presented in Figs. 9 and 10.All values and uncertainties are presented in Supplementary Information.
For the water-rich region of the curve, the behavior of the system was well-represented and predictions close to experimental data were obtained.For the methane-rich region, although some bigger deviations can be seen, the overall pressure-composition behavior was correctly represented, which can be considered satisfactory, especially considering the particular complexity of the system.Here is convenient to remember that the electronegativity difference between oxygen and hydrogen promotes intense polarization effects and hydrogen bonds, what leads to the formation of association structures.Moreover, these structures can have a variable number of water molecules.Hence, the correct representation of the water molecule and their association structures in the bulk is not trivial, especially with rigid potentials (Shvab 2014).The utilization of non-rigid potentials, though, increase simulation times.The challenge increases even further when far from ideal mixtures are simulated.Even a heavily parameterized equation of state like GERG-2008 presents some di culty in keeping the prediction accuracy seen before, for less challenging systems, like for CH 4 + CO 2 .
The underestimation of the water content in methane, though, is consistent with what was veri ed by Errington et al. (1998).The same was observed for water-ethane solubility.Bolton et al. (2009) evaluated the water solubility in decane and found the same trend for the classic Lorentz-Berthelot mixing rule and many of its modi cations, especially for temperatures under 450 K.Moreover, their simulations showed the same behavior is found for longer hydrocarbons, up to 300 carbons.
The comparison between the densities provided by equations of state and Monte Carlo simulations for the CH 4 + H 2 O system is given in Fig. 11, showing again an almost perfect agreement between simulations and equation of state.The exact values and uncertainties are presented in Supplementary Information.

Conclusion
Binary vapor-liquid equilibria of natural gas main components (CH 4 , CO 2 and H 2 O) were simulated through Monte Carlo simulations using the DSF method for the calculation of electrostatic interactions to check whether this less rigorous approach than the Ewald method can reliably predict natural gas systems' equilibria.For all studied systems, predictions comparable to experimental data were found.Furthermore, when confronting DSF simulations to other simulations with different approaches, including those with Ewald methods, DSF simulations showed comparable, or even superior, prediction accuracy.
It was also possible to see that the chemical complexity of the system in uenced directly in the prediction quality and indirectly in the simulation times.Thus, for CH 4 + CO 2 binary, high quality predictions were readily obtained, whereas for CO 2 + H 2 O more computational resources were required.Finally, the simulation times were increased multiple times for the last binary, CH 4 + H 2 O, and the obtained results were the less accurate among all studied systems, related with the complexity of the water charge distribution, its associating behavior and how it interacts with methane leading to very low mutual solubilities.
The results contribute to enable the DSF as a method that can aid the search for new solutions in the natural gas industry and to establish this approach as an alternative to make molecular simulations more popular and accessible, decreasing the computational demand and making homemade simulations for complex systems one step closer to reality.

Declaration of interests
The authors declare that they have no known competing nancial interests or personal relationships that could have appeared to in uence the work reported in this paper.

5 .
Arai et al. (1971) reported the experimental density behavior for this system at 253 and 273 K.All the data shown in Figs.4 and 5seen to be in qualitative and quantitative agreement with the values presented by Arai et al.

Figure 1 Example
Figure 1