Solubility of Carbon Dioxide, Hydrogen Sulfide, Methane, and Nitrogen in Monoethylene Glycol; Experiments and Molecular Simulation.

Knowledge on the solubility of gases, especially carbon dioxide (CO2), in monoethylene glycol (MEG) is relevant for a number of industrial applications such as separation processes and gas hydrate prevention. In this study, the solubility of CO2 in MEG was measured experimentally at temperatures of 333.15, 353.15, and 373.15 K. Experimental data were used to validate Monte Carlo (MC) simulations. Continuous fractional component MC simulations in the osmotic ensemble were performed to compute the solubility of CO2 in MEG at the same temperatures and at pressures up to 10 bar. MC simulations were also used to study the solubility of methane (CH4), hydrogen sulfide (H2S), and nitrogen (N2) in MEG at 373.15 K. Solubilities from experiments and simulations are in good agreement at low pressures, but deviations were observed at high pressures. Henry coefficients were also computed using MC simulations and compared to experimental values. The order of solubilities of the gases in MEG at 373.15 K was computed as H2S > CO2 > CH4 > N2. Force field modifications may be required to improve the prediction of solubilities of gases in MEG at high pressures and low temperatures.

Introduction

Monoethylene glycol (MEG) is a colorless, low-volatility, and stable liquid. MEG is fully miscible in water as well as in many organic liquids such as acetone and methanol.[1] In 2020, the global market size of MEG is valued at USD 20 billion.[2] MEG is widely used as an antifreeze agent, coolant, and heat-transfer agent and as a raw material for the manufacturing of polyester fibers.[3] In the oil-and-gas industry, MEG is widely used for the prevention of gas hydrate formation.[3,4] In the course of mitigating gas hydrate formation, MEG has been reported to absorb acid gases such as carbon dioxide (CO2) and hydrogen sulfide (H2S).[5] Because of the absorption capability, stability, and miscibility of MEG in many organic liquids, it is also considered for use in separation processes for acid gases.[5−8]

A number of MEG-based solvents, such as deep eutectic solvents, are considered for CO2 capture.[8−12] More recently, mixtures made from MEG, amines, and water are investigated for simultaneously preventing hydrate formation and removing H2S in offshore oil-and-gas applications.[13] To achieve these purposes, triethylene glycol (TEG)–amine–water mixtures were previously used. Replacing TEG with the less viscous MEG is expected to improve the absorption capability of glycols–amine–water solvents because absorption rates increase with lower viscosities.[14] To design and optimize processes in which MEG acts as a hydrate formation inhibitor or as an absorbent, knowledge of the phase equilibria of the system is essential.[7,14] To this purpose, a number of experimental measurements of binary mixtures of MEG and gases, that is, CO2 and H2S, have been performed.[6,7,15] For a review of experimental studies on the solubility of acid gases in MEG, the reader is referred to ref (7).

While traditionally phase equilibrium data are obtained from experimental measurements, such an approach is not always feasible, especially if high pressures and/or temperatures are required, and when dangerous gases, such as H2S, are involved. For this reason, theoretical approaches for computing the phase equilibria of mixtures of gases and liquids have been widely used.[16−19] Unlike classical thermodynamic models, molecular-based methods account for the strong molecular interactions present in associating liquids.[5] In the past few decades, molecular simulation has emerged as a powerful tool for using microscopic information of associating liquids to predict their macroscopic behavior.[20,21] In addition to providing thermodynamic and transport data, molecular simulation can also be used to investigate the microscopic structure of solvents and to understand absorption mechanisms.[22,23]

Monte Carlo (MC) simulations have been used to predict the solubility of gases in associating liquids.[22,24,25] MC has also been widely used to study the absorption of gases in solvents such as alcohols,[25−27] ionic liquids,[28,29] and deep eutectic solvents.[30] To the best of our knowledge, studies reporting MC simulations of the phase equilibria of small gas molecules and MEG are lacking. A possible reason for the absence of such studies is the fact that the simulation of dense liquids with strong intermolecular interactions, as in the case of MEG, is computationally demanding. MC simulations in open ensembles are often used to compute the solubility of solutes in liquids. In these ensembles, the solutes are added to or removed from the simulation box. For dense liquids and/or with the presence of strong interactions, such insertions can be challenging.[31,32] To enhance the efficiency of molecular transfers in MC simulations, Shi and Maginn[24,33] developed the continuous fractional component MC (CFCMC) simulation method. In this method, the system is expanded using a so-called fractional molecule with a coupling parameter λ which is used to vary the interactions between the fractional molecule and the surrounding molecules. In solubility calculations, a fractional molecule is used to gradually add/remove molecules to/from the solvent.[34] The presence of a fractional molecule does not affect the prediction of thermodynamic properties of the system.[32,35] For an in-depth discussion on the CFCMC method, the reader is referred to the recent review by Rahbari et al.[36]

A prerequisite for successful MC simulations of pure and multicomponent mixtures is the use of force fields that can adequately represent inter- and intramolecular interactions. Thus, another challenge of simulating gases in associating liquids, such as MEG, is the availability of force fields that provide accurate predictions of the desired properties. One of the most commonly used force fields for a large number of gases and liquids is the transferable potentials for phase equilibria (TraPPE).[37,38] The TraPPE force field has been successfully used for the prediction of thermodynamic and transport properties of gases and liquids.[25,27,28,39−41,41−44] Cardona et al.[45] used TraPPE and other classical force fields to compute the thermodynamic properties of pure MEG. The authors found that the united-atom version of TraPPE (TraPPE-UA)[38] is able to accurately predict the thermodynamic properties of pure MEG, such as the density, isothermal compressibility, and heat of vaporization.

The objective of this paper is to investigate the ability of MEG to absorb various gas molecules relevant to industrial applications. To study the solubility of CO2, calorimetric experiments are performed. MC simulations in the osmotic ensemble are carried out using the CFCMC method.[32,34] CFCMC calculations are compared to experimental measurements of CO2 in MEG. In addition to CO2, we study the absorption of other gases relevant to the oil-and-gas industry. MC simulations are used to predict the solubilities of CH4, N2, and H2S in MEG. The TraPPE-UA force field is used for MEG and gases studied in this work. In the case of the solubility of H2S, we test the performance of a force field developed by Kristóf and Liszi,[46] which was used in a number of solubility studies.[47,48] For the solubilities of CH4, N2, and H2S in MEG, the predictions of MC simulations are compared with experimental data from the literature and the performance of the force fields used is evaluated. For all systems, the Henry coefficients are computed using CFCMC simulations.

The paper is organized as follows: in Section , details related to the experimental setup are provided. In Section , the MC simulation methods used to compute the solubilities in MEG are explained. In Section , the experimental data and MC calculations of absorption isotherms of CO2 in MEG are provided. MC simulation results of the solubility of CH4, N2, and H2S in MEG are shown in Section . In Section , the main findings of this work are summarized.

Experimental Details

In this section, we describe the experimental setup used to measure the solubility of CO2 in MEG. In Table , the CAS registry number and mass fraction purity of MEG and CO2 are provided. The apparatus used to measure CO2 solubility in MEG consists of a CPA 122 calorimeter purchased from ChemiSens AB. This equipment is described in the previous works by Kim et al.[49,50] A schematic representation of this setup is shown in Figure . Here, we describe the experimental setup while referring to the number labeling (shown in the parentheses) of each device in Figure . Essentially, the setup consists of a large stirred cell reactor (1) of V ≈ 2 dm3 and two large cylinders (3) for storing CO2 of V ≈ 4.55 dm3. The reactor is vacuumed and partially filled with ca. 1.2 dm3 of the solvent at the start of each experiment. Vacuuming is performed three more times at room temperature so that only ethylene glycol vapor in equilibrium with the liquid is kept in the cell. Meanwhile, the CO2 cylinders are kept in a water bath (4) coupled with a temperature controller so that their conditions are also supervised. The temperature controller (2) is set to increase the temperature of the reactor up to a desired setpoint, and the stirrer is turned on at approximately 500 rpm. After the setpoint is reached and equilibrium is achieved for a minimum of 30 min, the first CO2 injection can be performed. A computer (5) is used to send a signal to the mass flow controller (MFC) (6) allowing its opening for a short amount of time. There are temperature (7) and pressure (8) measuring devices coupled to the cylinders, and another one of each coupled to the stirred cell (9, 10). With data gathered both before and after each injection and with knowledge of the volumes of the cylinders and of the stirred cell, one is able to calculate the amount of CO2 transferred to the gas phase and to the liquid phase of the reactor. For these calculations, the Peng–Robinson equation of state (PR-EOS) has been employed to correlate pressure–volume–temperature (PVT) to the number of moles (n). The PR-EOS requires the critical pressure Pc = 73.90 bar, the critical temperature Tc = 304.21 K, and the acentric factor ω = 0224, all of these values referring to pure CO2. No binary interaction parameters are used for these calculations. Additionally, the volumetric flow of CO2 passing through the MFC (6) is recorded by the computer (5) so that the amount of gas leaving the cylinders can be calculated either by PR-EOS or by simply reading the MFC data log. After having assessed that both methods return roughly the same values, all solubility results shown in this work have been obtained by reading the values recorded by the MFC.

Table 1

CAS Registry Number and Mass Fraction Purity of Components Used in the Experimental Measurements

component	CAS	supplier	mass fraction
MEG	107-21-1	Sigma-Aldrich	0.998
carbon dioxide	124-38-9	AGA	0.99999

Figure 1

Setup of the CPA 122 calorimeter used to measure the solubility of CO2 in MEG. MEG is filled in a stirred cell reactor (1) linked to a temperature controller (2). CO2 is stored in two cylinders (3) that are placed in a water path (4). A computer (5) sends signals to the MFC (6) to inject CO2 in the cell reactor (1). Temperature and pressure data of CO2 are measured using devices 7 and 8, respectively. For MEG, devices 9 and 10 are used. The power transferred from the reactor (1) to the temperature controller is measured using device 11.

The calorimeter CPA 122 is also able to record heat of absorption data from measuring (11) of the power transferred from the reactor to the temperature controller. Treatment of heat of absorption data was not performed because they are not relevant to this study. After each injection and after equilibrium is attained in the CO2 cylinders, these cylinders have to be refilled with CO2 from the wall (12) through opening of the valve (13). After achieving equilibrium in both cylinders and in the stirred cell, the next injection is performed and the process is repeated until the reactor approaches the maximum operational pressure of 6 bar.

The calorimeter CPA 122 does not directly deliver CO2 partial pressure data but total pressure data. After each injection, it is expected that the increase in pressure measured in the reactor is solely due to an extra amount of CO2 which is transferred to the gas phase and not due to any volatilization of MEG brought by the CO2 injection. In other words, treatment of the solubility data requires a constant partial pressure of MEG throughout the whole experiment. There are no grounds to suspect that this assumption is unreasonable. Therefore, the partial pressure of CO2 can be back-calculated by subtracting the partial pressure of MEG at the start of the experiment from the total pressure obtained after each injection.

The expanded uncertainties with a 0.95 level of confidence for the CO2 solubility experiments have been evaluated as being approximately U(α) = 0.1, implying that CO2 loadings are reported within a confidence margin of roughly 10%. A detailed procedure for the calculation of expanded uncertainties is provided in the Supporting Information of this study.

For these experiments, MEG anhydrous with a purity of 99.8% (CAS 107-21-1) has been supplied by Sigma-Aldrich and used in the procedures without further purification, while carbon dioxide of a 99.999% purity has been supplied by AGA.

MC Simulation Details

Force Fields

Classical force fields were used to describe the interactions of the molecules studied in this work. For MEG, all interaction potentials and parameters follow from the TraPPE-UA force field.[38] The TraPPE force field adequately predicts densities and vapor–liquid equilibria (VLE) of many species such as normal alkanes,[37] branched alkanes,[51] glycols, and ketones.[38] To accurately represent the molecular structure of MEG, Stubbs et al.[38] added an additional repulsive term (r–12) for interactions between a hydroxyl hydrogen and an oxygen atom separated by four bonds. In our study, the TraPPE-UA force field was used also to represent CO2, CH4, H2S, and N2 as rigid objects. We also tested another four-site model presented by Kristóf and Liszi[46] for H2S, referred to here as H2S-KL. The main differences between the two force fields are the nonbonded Lennard-Jones (LJ) parameters and the atomic charges. Table lists the components simulated in this study and the force field used for each component. All force field parameters are listed in the Supporting Information.

Table 2

Chemical Formulas and Force Fields of the Components Simulated in This Work

component	chemical formula	force field
MEG	HO(CH2)2OH	TraPPE-UA[38]
carbon dioxide	CO2	TraPPE-UA[70]
methane	CH4	TraPPE-UA[37]
hydrogen sulfide	H2S	TraPPE-UA[69]
hydrogen sulfide	H2S	Kristóf and Liszi[46]
nitrogen	N2	TraPPE-UA[70]

In this work, two types of intermolecular interactions are computed: LJ and Coulombic interactions. LJ interactions were truncated at 12 Å, and the uncertainty due to truncation is handled by applying analytic tail corrections.[20,21] The Lorentz–Berthelot mixing rules were used for LJ interactions between dissimilar interaction sites.[20,21] The Ewald summation method was applied to handle electrostatic interactions with a relative precision of 10–6. The real-space part of the electrostatic interactions was truncated at 12 Å. All simulations were carried out in the osmotic ensemble (see Section ). The PC-SAFT equation of state was used to compute the fugacity of the solutes at the desired temperatures and pressures.[52,53]

CFCMC Method

The osmotic ensemble[21] is used to compute the solubility of small solute molecules in nonvolatile solvents. In this open ensemble, the following parameters are fixed: the temperature (T), the hydrostatic pressure (P), the number of molecules of MEG (NMEG), and the fugacity of the solute (f). The number of molecules of the solute (Ns) and the volume of the system (V) are varied to achieve equilibrium. The hydrostatic pressure P inside the simulation box corresponds to the imposed fugacity f of the reservoir. An essential part of the calculations is the insertion and/or deletion of solute molecules in the simulation box. When studying dense solvents, as in the case of MEG, molecule insertions can be challenging.[31,32] To improve the probability of accepting insertion/deletion moves, the CFCMC method was used.[24,32,36,54,55] The osmotic ensemble was expanded using a so-called fractional molecule.[24] As opposed to a whole molecule, the strength of interactions of a fractional molecule is varied using a coupling parameter λ. When λ = 0, the fractional molecule acts as an ideal gas molecule and does not interact with the surrounding molecules. When λ = 1, the fractional molecule fully interacts with the surrounding molecules. By varying the strength of interactions of the fractional molecule with the surrounding molecules, whole molecules can be gradually added or removed. Besides the conventional MC thermalization trial moves,[20,21] trial moves attempting to change λ are required. Shi and Maginn[24] derived Metropolis-like acceptance rules for changing the values of λ in the osmotic ensemble. For each solute type, a fractional molecule is used to insert/delete molecules in the simulation box. When λ drops below 0, the fractional molecule is removed and a randomly selected whole molecule is transformed into a fractional molecule. When λ is larger than 1, the fractional molecule is transformed into a whole molecule and a new fractional molecule is inserted.[24,33] For the solvent, a fractional molecule is used to improve sampling, while keeping the total number of molecules of the solvent fixed. For a fractional molecule of the solvent type, λ trial moves involve changing the values of λ, random reinsertions of the fractional molecule, and identity swaps between a fractional molecule and a whole molecule.[35] In CFCMC simulations, the system is biased to ensure a flat distribution of the observed probability of λ. From the probability distribution of λ, the excess chemical potential and hence the Henry coefficient are obtained. For more details, the reader is referred to refs.[24,32,34,35]

Scaling of the Intermolecular Interactions

For each fractional molecule, a weight function W(λ) is constructed to achieve a flat probability distribution of λ and ensure smooth transitions between λ = 0 and λ = 1.[32] Essentially, at a certain λ, the value of W(λ) counteracts the value of ⟨∂U/∂λ⟩, which is the average potential energy change with λ. Fluctuations in the value of ∂U/∂λ can be large, which can hinder the efficient sampling of the λ-space. As a result, a pathway that minimizes the variance of ∂U/∂λ has to be chosen.[56] Changes in the values of λ depend on how intermolecular interactions are scaled when λ is varied[56] (intramolecular interactions are not scaled). Electrostatic interactions are scaled by using a scaling parameter λel that is a function of λ.[34,57] For nonbonded LJ interactions, the following soft-core scaling potential is used[56,58,59]where i and j are the interaction sites, ϵ and σ are the LJ parameters, and r is the distance between i and j. The parameters a, b, c, and α are adjusted to achieve an efficient sampling of the λ-space. For systems composed of MEG molecules, a number of scaling potentials were tested. Figure shows the values of λ of an MEG fractional as a function of the number of MC cycles for three scaling potentials. The commonly used 1-2-6 potential results in poor sampling of the λ-space. Figure a shows that at certain periods, the values of λ are confined to a limited range. Changing the parameter b from 2 to 1 improves the sampling as demonstrated in Figure b. Figure c shows that the 1-1-48 potential with α = 0.0025, which was recommended earlier by Pham et al.,[56] results in the best sampling.

Figure 2

Values of λ vs the number of MC cycles of an MEG fractional molecule in NPT MC simulations. To scale the interaction of the fractional molecules, the following scaling potentials are used (see eq ): (a) 1-2-6, (b) 1-1-6, and (c) 1-1-48.

Simulation Details

Molecular simulations were performed using the recently developed open-source software package Brick-CFCMC.[34] The density of pure MEG was computed in the NPT ensemble at 1 bar and at temperatures of 333.15 and 353.15 K. The solubility of CO2 in MEG was computed at three temperatures, T = 333.15, 353.15, and 373.15 K. The solubilities of CH4, N2, and H2S in MEG were computed at T = 373.15 K. For all gases, solubilities were computed at pressures ranging from 1 to 10 bar, but in the case of N2, pressures up to 100 bar were considered because N2 has very low solubilities in MEG at low pressures. Simulation boxes were set up with 250–350 MEG molecules depending on the number of solute molecules absorbed. Two fractional molecules were used: a fractional molecule to insert/remove solute molecules into the simulation box and a MEG fractional molecule. Adding an MEG fractional to the simulation box has the following advantages: (1) it enhances sampling as the fractional molecule can be used for random reinsertions (low λ) and identity changes (high λ) and (2) the excess chemical potential of MEG is automatically computed, from which the saturated vapor pressure of MEG can be estimated. The following MC trial moves were used: translations, rotations, and volume change trial moves. MC trial moves that attempt to change the values of λ were used for both fractional molecules. Simulations in the osmotic ensemble were carried out with the following probabilities for selecting trial moves: 25% translations, 25% rotations, 32% intramolecular moves, 1% volume changes, and 17% λ trial moves (divided equally between the solute and MEG fractional molecules). A minimum of 1 × 106 cycles were carried out for equilibration. At each MC cycle, the number of the trial moves performed equals the number of molecules of the system.

During equilibration, an iterative scheme was used to obtain a weight function W(λ) that results in a flat probability distribution of λ. For production runs, 1 × 106 cycles were carried out. To minimize the statistical error of the computed averages, a number of independent production simulations were performed at a specified T and P. The number of simulations performed was selected such that the uncertainty is less than 5%. For each system simulated in this work, at least 25 independent production runs were carried out. Error bars were computed by dividing these runs into five groups and calculating the standard deviation.[21] For each system, these error bars are reported in the Supporting Information (Tables S14–S19).

Results and Discussion

Solubility of CO2 in MEG from Experiments and MC Simulations

Densities of pure MEG were computed using MC simulations in the NPT ensemble at P = 1 bar and at T = 333.15 K and T = 353.15 K. In Table , densities (reported in units of kg/m3) from simulations are compared to experimental values from the work of Skylogianni et al.[13]Table shows that when using the TraPPE-UA force field, simulations underpredict the densities of MEG. The differences between experiments and simulations are around 5%. Simulating a solvent with an underestimated density may result in higher absorption capacity. Deviations between experiments and simulations will be discussed in detail later in this section.

Table 3

Properties of Pure MEG at Different Temperatures from Experiments and MC Simulationsa

T [K]	ρexpL [kg/m3]	ρsimL [kg/m3]	μex/kB [K]	ρsimV [kg/m3]	Psimsat [bar]	PNISTsat [bar]
333.15	1085	1029.1 ± 20	–4125 ± 41	4.2 × 10–6	0.022	0.021
353.15	1070.1	1016.3 ± 0.5	–3948 ± 25	1.4 × 10–5	0.0068	0.0067
373.15		1003.3 ± 0.5	–3749 ± 27	4.3 × 10–5	0.0019	0.0018

Density of pure MEG in kg/m3 obtained from experiments[13] and MC simulations at P = 1 bar. The vapor densities and saturated vapor pressures of MEG are computed from μex/[kBT] (see Section ). Experimental vapor pressures are obtained from the NIST database.[60]

As a fractional MEG molecule is present in the simulation, we can calculate the excess chemical potential of MEG μex from the probability distribution of its λ parameter. The chemical potential of MEG in the liquid phase equals[35]where ρL is the number density of MEG and μo is the ideal gas chemical potential, which only depends on the temperature.[21]Equation also applies to MEG in the gas phase. At equilibrium, the chemical potentials in the liquid phase and gas phase are equal. If we assume an ideal gas phase, then μex in the gas phase equals 0 and ρV = Psat/kBT. From this, the saturated vapor pressure Psat of MEG can be estimated (eq. S7 in the Supporting Information). In Table , we report the excess chemical potential, vapor densities, and saturated vapor pressures of pure MEG. The vapor pressures computed from MC simulations are compared to experimental values obtained from the NIST database.[60]Table shows that both methods are in good agreement. The pressures reported in Table can be considered very small, which validates the assumption made by the experimental method regarding the nonvolatility of MEG.

Figure shows absorption isotherms of CO2 in MEG from experiments and MC simulations in the osmotic ensemble at the temperatures of T = 333.15 K, T = 353.15 K, and T = 373.15 K. In Figure , a typical MC simulation snapshot of MEG and CO2 molecules is shown. Figure shows that CO2 molecules are dispersed in MEG and not clustered. In Figure , the ratio of the number of moles of CO2 to the number of moles of MEG (i.e., loading) is plotted as a function of pressure. Solubilities from MC simulations were found to qualitatively agree with experimental measurements. Both experiments and simulations show that for all temperatures studied, the loading is almost linear in this pressure range. Additionally, both methods report that the absorption of CO2 in MEG decreases at higher temperatures. Figure shows that the quantitative agreement between MC simulations and experiments varies with temperature and pressure. At very low pressures (i.e., <2 bar), loadings of CO2 obtained from MC simulations agree well with experiments. At higher pressures, MC simulations overpredict the absorption of CO2 when compared to experiments. The deviation between simulations and experiments systematically increases with pressure and decreases with temperature. The inherent uncertainties of the experimental loadings of CO2 in MEG are shown with vertical error bars. Uncertainties of experimental values are calculated using a methodology described in the work of Wanderley and co-workers.[61] The inherent uncertainties of the total pressure are delimited by the sensitivity of the pressure transducer employed for those measurements, which is ±0.015 bar. Conversely, this implies that the uncertainty of the estimated partial pressures of CO2 is ±0.021 bar because of error propagation. One can see in Figure that these are very small uncertainties when considering the span of pressures measured in the series of experiments. In the Supporting Information, solubilities from experiments and MC simulations are provided in a tabulated from along with their uncertainties.

Figure 3

Absorption isotherms of CO2 in MEG at different temperatures. Closed symbols are solubilities using MC simulations (details in Section ), and open symbols are experimental results of this work. The raw data and the corresponding uncertainties are provided in Tables S7–S9 and S14–S16 of the Supporting Information.

Figure 4

(a) Typical snapshot of a simulation of MEG in which CO2 is absorbed in the osmotic ensemble (T = 333 K, P = 8 bar, NMEG = 220 molecules, and NCO = 3 molecules) in a simulation box with the dimensions of 28 × 28 × 28 Å. (b) Same snapshot as in (a) while showing only CO2 molecules. Clearly, the CO2 molecules are not clustered. Figures were produced using the software package iRASPA.[71]

In Figure , solubilities measured in this work at T = 373.15 K are compared to solubilities from other experimental studies. The measurements in this work were found to match the data from Jou et al.[15] at low pressures. The experimental method of Jou et al.[15] differs from our experimental method. The main difference is that Jou et al.[15] evaluated CO2 concentrations using gas chromatography and acid–base titration, while we apply a mass balance approach. At higher pressures, CO2 solubilities from other experimental works slightly differ from the results of this work. Figure also shows that loadings computed using MC simulations agree the most with our experimental results. In the Supporting Information (Figure S2), we compare our simulations and experimental results with other experimental data from the literature.[14,62,63]

Figure 5

Absorption isotherm of CO2 in MEG at T = 373.15 K. Closed symbols are solubilities using MC simulations (details in Section ), and open symbols are experimental results. The raw data and the corresponding uncertainties are provided in Tables S9, S10, and S16 of the Supporting Information.

Besides absorption isotherms, it is also possible to describe the solubility of gases in solvents through Henry coefficients. The Henry coefficient of solute 2 in solvent 1 is defined as[64]where P2 and x2 are the partial pressure and mole fraction of solute 2, respectively, and f2 is its fugacity. With these experimental values, the Henry coefficient is defined as the partial pressure of CO2 in bar divided by the CO2 molar fraction in MEG. In MC simulations, Henry coefficients H21 are computed from the excess chemical potential of the solute μ2ex[65]

In CFCMC, μ2ex is computed from sampling the probability distribution of λ.[32,35] In Table , Henry coefficients of CO2 in MEG HCO computed using MC simulations are reported at different temperatures and compared to Henry coefficients from experiments. Both methods demonstrate that the value of HCO increases with temperature. The maximum difference between experimental and computed Henry coefficients is 30%. The difference consistently decreases with increasing temperature to reach 13% at T = 373.15 K. Predictions from MC simulations are satisfactory considering that the force fields and the mixing rules used for MEG and CO2 were not modified. The Henry coefficients reported in Table indicate that pure MEG would not be a good absorbent for CO2. In a study by Ramdin et al.,[66] Henry coefficients at T = 333 K of CO2 in selexol and the ionic liquid [bmim][TF2N] were reported to be 68 bar and 66 bar, respectively. At the same temperature, the experimental Henry coefficient of CO2 in MEG is 634 bar (Table ).

Table 4

Henry Coefficient of CO2 in MEG at Different Temperatures Obtained from the Experiments (This Work) and Molecular Simulations Performed in This Work

	HCO2,MEG, bar (mol CO2/mol EG)
T [K]	experimental	MC simulations
333.15	634 ± 30	445 ± 20
353.15	736 ± 36	576 ± 15
373.15	843 ± 40	730 ± 21

From the knowledge of solubility of CO2 at different temperatures, the heat of absorption q of CO2 in MEG can be calculated usingwhere R is the gas constant and Po is a reference pressure to make the argument of the logarithm dimensionless. Using solubilities of CO2 in MEG from MC simulations of this work, q was found to be equal to −12.8 kJ/mol, indicating that the absorption of CO2 is an exothermic process. This value is in good agreement with experimental findings. Recently, the heat of absorption was measured using calorimetric experiments in a study by Wanderley et al.[67] and was reported to be −14 kJ/mol at 343.15 K.

The differences between theoretical and experimental solubilities can be attributed to the force field used to describe MEG. From our simulations, it is observed that the TraPPE-UA force field underpredicts the density of pure MEG (see Table ). Lower MEG densities can potentially lead to higher absorption capacities of solutes. Moreover, the force field parameters of TraPPE-UA[38] were obtained using the VLE experimental data of MEG at high temperatures (>400 K), and as a result, inaccuracies at lower temperatures can be expected as we move outside the fitting range of the TraPPE force field.

Figure shows that deviations between experiments and simulations are larger at lower temperatures. While deviations can be reduced by optimizing the force field parameters of MEG, force field parameters of the solute have to be considered as well. For CO2, TraPPE force field parameters are obtained using pure component data and not data of multicomponent systems. Predictions of MC simulations can be improved by revising force field parameterization or considering different force field combinations. Alternatively, one might consider changes to the combination rules used. Here, the Lorentz–Berthelot rules are used to compute interactions between dissimilar sites. To improve the predictions of MC simulations, other combination rules can be considered and/or adjustable parameter(s) can be added to fine-tune solute–solvent interactions.

Solubility of CH4, H2S, and N2 in MEG from MC Simulations

MC simulations were used to compute the solubility of other pure gases in MEG at 373.15 K. MC simulation results were compared to experimental data from the literature.

In Figure , the absorption isotherm of CH4 in MEG is shown for T = 373.15 K and pressures ranging from 1 to 10 bar. At this pressure range, low loadings of CH4 are obtained from MC simulations. To validate computational results, experimental solubilities[68] at higher pressures are shown in Figure . At P = 17.9 bar, MC simulations overpredict the solubility of CH4 in MEG by 25%. As discussed earlier in Section , higher absorption of solutes is due to the underestimated densities of MEG when using the TraPPE-UA force field.

Figure 6

Absorption isotherm of CH4 in MEG at T = 373.15 K. Closed symbols are solubilities from MC simulations in the osmotic ensemble (details in Section ), and open symbols are experimental data from ref (68). The raw data and the corresponding uncertainties are provided in Tables S11 and S17 of the Supporting Information.

In Figure , solubilities of H2S in MEG at T = 373.15 K from MC simulations using two different H2S force fields are compared to experimental solubilities from ref (15). A reasonable agreement between MC simulations and experiments is obtained for the two force fields, but larger deviations increase at high pressures. The H2S-TraPPE force field underpredicts loadings of H2S, while the H2S-KL force field overpredicts loadings under the studied conditions. At atmospheric pressure, solubilities computed using H2S-KL were found to be closer to the experimental value reported by Jou et al.,[15] compared to the solubility computed using H2S-TraPPE.

Figure 7

Absorption isotherm of H2S in MEG at T = 373.15 K. Closed symbols are solubilities from MC simulations (details in Section ) using two force fields: H2S-TraPPE and H2S-KL.[46] Open symbols are experimental data from ref (9). The raw data and the corresponding uncertainties are provided in Tables S13 and S18 of the Supporting Information.

In Figure , the absorption isotherm of N2 in MEG at T = 373.15 K is shown. The loading of N2 is computed using MC simulations and is compared to experimental data from Zheng et al.[6] Because the absorption of N2 in MEG is negligible at atmospheric pressures, simulations were performed at pressures ranging from 10 to 100 bar. From Figure , it can be seen that the computed loadings deviate considerably from experimental values and that the deviations increase systematically with pressure. As discussed earlier, differences between MC simulations and experimental data at high pressures can be improved by modifying the used force fields or fine-tuning solute–solvent interactions.

Figure 8

Absorption isotherm of N2 in MEG at T = 373.15 K. Closed symbols are solubilities from MC simulations (details in Section ), and open symbols are experimental data from ref (6). The raw data and the corresponding uncertainties are provided in Tables S12 and S19 of the Supporting Information.

In Table , Henry coefficients computed using MC simulations of CH4, H2S, and N2 are listed. Experimental Henry coefficients of CH4 and H2S are also shown. The average differences between experimental and computational values are around 25%. From Henry coefficients at T = 373.15 K, the following order of solubility is exhibited: H2S > CO2 > CH4 > N2. From Henry coefficients, the ideal selectivity of the desired component i in the undesired component j can be expressed as

Table 5

Henry Coefficients of Various Solutes in MEG Obtained from Experiments and MC Simulationsa

	Hi,j, bar (mol i/mol EG)
solute	experimental	MC simulations
CH4	5673	4504 ± 165
H2S-TraPPE	227.4	302 ± 14
H2S-KL	227.4	173 ± 5
N2		10,815 ± 248

Experimental values are taken from refs (6) and[68] for CH4 and H2S, respectively.

In Table , ideal selectivities for the separation of CO2 using MEG are provided, which are computed using Henry coefficients from experiments and MC simulations. The results show that CO2 is more selectively absorbed in MEG than CO4 and N2. However, this is not true for the separation of CO2 using MEG in the presence of H2S.

Table 6

Ideal Selectivities of CO2, CH4, H2S, and N2 in MEG Computed Using Henry Coefficients from Experiments and MC Simulations at T = 373.15 Ka

	Si/j
separation	experimental	MC simulations
CO2/CH4	6.73	6.17
CO2/H2S-TraPPE	0.27	0.41
CO2/H2S-KL	0.27	0.24
CO2/N2		14.81

Experimental values are taken from refs (6) and[68] for CH4 and H2S, respectively.

Conclusions

In this work, the solubility of CO2 in MEG was studied both experimentally and computationally. MC simulations in the osmotic ensemble were performed to predict the absorption of CO2 as well as the absorption of CH4, H2S, and N2. The CFCMC method was used to facilitate the insertion/deletion of particles into the solvent. TraPPE force fields were used to model all species. For H2S, two force fields were compared: H2S-TraPPE[69] and H2S-KL.[46] The solubility of CO2 in MEG was measured at the following temperatures: 333.15, 353.15, and 373.15 K. From experiments and MC simulations, CO2 was found to be better absorbed at lower temperatures. At T = 373.15 K, CO2 and H2S were found to have higher solubilities in MEG than CH4 and N2. Solubilities predicted by MC simulations are in reasonable agreement with experimental data. For all the solutes studied in this work, deviations between MC simulations and experiments were found to increase with pressure. For the solubility of H2S, predictions from the H2S-KL force field were closer to experimental measurements than those from H2S-TraPPE. Other than absorption isotherms, Henry coefficients were also computed. The order of solubilities of the gases in MEG at 373.15 K was found to be as follows: H2S > CO2 > CH4 > N2. The average difference between Henry coefficients from experiments and Henry coefficients from MC simulations is around 20%. These results can be regarded satisfactory, considering that force fields from the literature were directly used without fitting binary interaction parameters. To improve predictions at high pressures, force field adjustments are required. For the solubility of CO2 in MEG, the experimental data provided in this work may be used to generate new force field parameters.