Pre-Screening of Ionic Liquids as Gas Hydrate Inhibitor via Application of COSMO-RS for Methane Hydrate

Abstract

Ionic liquids (ILs) due to their potential dual functionality to shift hydrate equilibrium curve and retard hydrate nucleation are considered as a very promising gas hydrate inhibitor. However, experimental testing alone is insufficient to examine all potential ILs combinations due to a high number of cation and anion to form ILs. In this context, four fundamental properties of IL-hydrate system, namely, sigma profile, hydrogen bonding energies, activity coefficient, and solubility, were stimulated through conductor-like screening model for real solvent (COSMO-RS). ILs were then analyzed to determine if they can be correlated with IL inhibition ability. Among them, sigma profile and hydrogen bonding energies, which later upgraded to total interaction energies, exhibit a significant relationship with IL inhibition ability. Total interaction energies of ions, on the other hand, have successfully been applied to develop a model. The model can predict the thermodynamic inhibition ability in terms of average temperature depression. The correlation was further validated with experimental values from literature with an average error of 20.49%. Finally, using sigma profile graph and developed correlation, the inhibition ability of 20 ammonium-based ILs (AILs) have been predicted. Tetramethylammonium hydroxide (TMA-OH), due to its short alkyl chain length cation and highly electronegative anion, has shown the most promising inhibition ability among the considered system.

1. Introduction

Gas hydrates are icelike crystalline solid compounds that could form in the presence of water and gas under favorable thermodynamic temperature-pressure condition [1]. At low-temperature and high-pressure conditions, water molecules (host) will surround the gas molecules (guest) and encapsulate the gas in a hydrogen-bonded solid lattice [2]. Depending on the gases trapped, different structures of gas hydrates can be formed. The structure I hydrate trapped methane (CH₄), ethane (C₂H₆), and carbon dioxide (CO₂) gases. Structure II usually forms for propane (C₃H₈) gas, while a mixture of CH₄ and butane (C₄H₁₀) and other hydrocarbons can be captured by structure H hydrates [3].

In recent decades, hydrates have received plenty of attention, because of its potential to capture and store gas [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21]. Also, it is discovered that gas hydrates located in subsea as well as permafrost region are a potential source of energy too [15]. However, the formation of natural gas hydrates in oil and gas pipeline is never applauded [4, 15, 18, 22]. This is because hydrate formation in pipelines has resulted in blockage and affected flow assurance of natural gas [23]. In spite of the economic losses caused by the blockage, ecological disasters could occur in severe cases too [24]. To prevent hydrate formation, several methods including isobaric thermal heating, water removal, depressurization, and chemical inhibitor injection [25] have been implemented. The three former methods, however, are not feasible and costly. As a result, chemical inhibitors have been researched and developed a lot in recent years to control the growth of hydrates.

There are generally three types of inhibitors, which are thermodynamic hydrate inhibitor (THI), kinetic hydrate inhibitor (KHI), and anti-agglomerates (AA). THI prevents the formation of the hydrate by shifting the thermodynamic equilibrium curve of gas hydrate to a lower temperature and higher pressure [25]. KHI, on the other hand, does not inhibit hydration formation, but it slows down their nucleation and growth of hydrate. It works on the principle of lengthening the formation time of hydrate to be longer than the residence time of the gas in pipelines [26]. Finally, AA, also a low-dosage inhibitor, allows the formation of hydrate but, through perturbation of water molecules, prevents the hydrate molecules from accumulating and growing larger [27].

Some common THI inhibitors include methanol and sodium chloride. To be effective, THI normally needs to be injected in a high concentration of around 10–50 wt% [28], which leads to high operational cost. Furthermore, sodium chloride corrodes oil and gas pipelines [29]. While KHI inhibitors were able to work effectively at a lower dosage (<1 wt%), Kelland reported that as exploration operation goes into the deeper sea, KHI still has to work together with THI to effectively inhibit hydrate formation [27]. These limitations signify that existing chemical inhibitors are still not performing well, and there is a strong need to develop more effective inhibitor [29, 30].

This leads the oil and gas industry toward ILs which was initially introduced as inhibitors by Chen et al. [31] in 2008, as the team discussed the effect of 1-butyl-3-methylimidazolium tetrafluoroborate in inhibiting CO₂ hydrate formation. A year later, Xiao and Adidharma [29] suggested the dual function of ILs inhibitors. The results showed that IL is not only able to shift the hydrate thermodynamic equilibrium curve, but it also retards the formation of the hydrate. Since then, numerous experimental works have been carried out to study the effect of ILs in inhibiting gas hydrates formation, mainly using imidazolium- and pyridinium-based ILs [2, 3, 25, 32]. The targeted ILs of this context are ammonium-based ILs (AILs), which are cheaper and easier to synthesis, but not being studied intensively. Therefore, due to cost economics and more environmentally friendly, AILs are chosen to be studied in this work.

To date, all the testing work of ILs effectiveness is done using an experimental method, which is by measuring the average depression temperature for thermodynamic hydrate inhibitors and by measuring induction time for kinetic hydrate inhibitors. There are generally no other methods available to validate the experimental work or to pre-screen ILs in a shorter time. Due to this reason, it is very desirable if a theoretical method to predict ILs effectiveness as hydrate inhibitors could be established just by analyzing their fundamental properties. And to obtain these fundamental properties, COSMO-RS, a thermodynamic properties predictive tool, is the best option available in the market.

For this purpose, COSMO-RS, which can estimate the fundamental properties of ILs system, has been selected. COSMO-RS is a novel method to predict the thermodynamic properties of ILs based on quantum chemistry model [33]. COSMO-RS first calculates the charge density of individual molecules based on the structure of each molecule [34]. The charge density will then be distributed onto the entire molecule surface. This distribution will then be described by a one-dimensional probability density [35], or more famously known as sigma profile, P(σ). Lastly, from the charge density, chemical potential, μ, will be calculated, and it will act as the basis for all other calculations to predict thermodynamic properties such as Henry’s law constant and activity coefficient [36]. The calculated properties will then try to be correlated to IL inhibition ability to develop a prediction model that could predict the inhibition ability of ILs.

Throughout the years, COSMO-RS model has been successfully applied in numerous works to predict the thermodynamic properties of systems containing ILs, such as liquid-liquid equilibrium [37, 38] and activity coefficient [34, 39]. Therefore, this has prompted a lot of screening efforts of ILs through COSMO-RS for different purposes such as determining extraction solvent and improving separating process [37, 40, 41, 42]. Grabda et al. [43], for example, has used COSMO-RS to carry out a screening process for ILs that is used as an extraction solvent for neodymium chloride and dysprosium chloride. Kurnia and Mutalib [44], on the other hand, had screened imidazolium-based ILs for the separation process of benzene from n-hexane through COSMO-RS. Other than screening work, comparison and validation work have been conducted too. Calvar et al. [37], for instance, have compared COSMO-RS prediction of LLE values of ILs with their experimental data and found out that the result is satisfactory. In 2007, Palomar et al. [45] reinforced the applicability of COSMO-RS in predicting density and molar volume of imidazolium-based IL when their predicted values laid close to the experimental data.

To support the application of COSMO-RS in this work, it is found out that many other applications involving ammonium-based and bionic ILs have already been conducted through COSMO-RS [43]. In 2010, Sumon and Henni [46] performed a COSMO-RS study on the properties of ILs for CO₂ capture. In this study, 12 ammonium-based cations such as tetramethylammonium (TMA), tetraethylammonium, and tetrabutylammonium (TBS) cations are used to derive ammonium-based ILs to be studied. In 2014, Grabda et al. [43] studied the effectiveness of 4400 ILs for NdCl₃ and DyCl₃ extraction. Among the many cations used are tetra-n-butylammonium, tetraethylammonium, tetramethylammonium, etc. Dodecyl-dimethyl-3-sulfopropylammonium cation, which is a type of ammonium-based cation, was concluded as the best performing cation in decreasing the chemical potential of NdCl₃ and DyCl₃, thus increasing their solubility and easing the extraction process. In the same year, Pilli et al. [47] screened out the best ILs to extract phthalic acid from aqueous solution using COSMO-RS. Although ammonium-based cation ILs in this simulation do not give the highest selectivity, they, however, have the highest activity coefficient. Next, through COSMO-RS, Machanová et al. [48] also obtained well-predicted values of excess molar volumes and excess enthalpy for N-alkyl-triethylammonium-based ILs.

As it observed from literature, screening of ILs for gas hydrate inhibition through COSMO-RS is a relatively new and fresh concept, yet, based on the successfulness of previous works [6, 18, 21, 22, 37] in predicting thermodynamic properties which provide the way for this current work.

2. Methodology

The research methodology comprises several activities described in even detailed and specific manner.

2.1 Extracting experimental IL inhibition ability

As a relatively new study, it is very important to gain acknowledgment and recognition from peers. Hence, as mentioned earlier, the experimental value of IL inhibition ability will be obtained from several past studies that are highly recognized. For instance, paper from Xiao et al. is chosen as it is the pioneer of IL inhibitor research. The full list of papers that were chosen for development or correlation and later for validation work is shown in Table 1.

Table 1.

Chosen papers for experimental values for this work.

As observed from table, experimental values from four papers will be collected. All of them studied hydrate formation in the presence of methane gas for the thermodynamic hydrate inhibitor. In all these papers, the effectiveness of an IL as THI was reported in the form of IL-hydrate equilibrium curve. Generally, a larger temperature depression signifies that the IL is good in inhibiting and shifting the equilibrium curve. However, since IL-hydrate equilibrium curve is not quantifiable and thus is not possible to develop correlation, average temperature depression will be used to represent IL inhibition ability in our work. This average temperature depression value can be calculated through the following equation [9, 13, 30]:

$$\stackrel{\mathrm{̶}}{T}=\frac{\sum \mathit{\Delta }\mathit{T}}{n}=\frac{{\sum }_{i=1}^n\left(T_{0,\mathrm{pi}}-T_{1,\mathrm{pi}}\right)}{n}$$ E1

where T_{0, pi} is the dissociation temperature of methane in a blank sample without IL and T_{1, pi} is the dissociation temperature of methane in a sample with IL inhibitor. The values of both dissociation temperatures should be obtained from the same p_i, and n refers to the number of pressure point considered. For example, Figure 1 shows the IL-hydrate equilibrium curve from Keshavarz et al. [49] for blank hydrate system (without IL) and hydrate system with 1-butyl-3-methylimidazolium tetrafluoroborate (BMIM-BF₄). Now, it is seen that with IL that acts as an inhibitor, the region of hydrate formation has reduced. It is also clear that the favorable pressure for hydrate to form has increased and the favorable temperature has reduced. This, in turn, made it hard for hydrate to form. Now to calculate average temperature depression, for instance, at 4 MPa, T_{0, pi} is equal to the temperature of blank hydrate without IL; the temperature would be around 277.5 K. On the other hand, T_{1, pi} that refers to the temperature of IL-hydrate system will be around 277 K. The difference between these two values is then the temperature depression. Several temperature depression values will be collected at different pressure points along the curve. Lastly, the average of these values will become the average temperature depression value.

Figure 1.

Hydrate-IL equilibrium curve from the work of Keshavarz et al.

2.2 Simulation of fundamental properties value in COSMO-RS

After obtaining the data of IL inhibition ability, now it is the time to collect another set of data, which is the fundamental property value of IL-hydrate system. Here, COSMO-RS software will be used to carry out the simulation. In COSMO-RS, all calculation works are performed based on density functional theory (DFT), utilizing the triple-zeta valence polarized (TZVP) basis set [50]. Figure 2 shows the entire computational method of COSMO-RS.

Figure 2.

Flowchart of predicting thermodynamic properties through COSMO-RS.

As regards Figure 2, COSMO-RS first requires the input of molecular structure [51]. After this, the charge density of a segment on each molecule surface will be calculated in a virtual conductor. The distribution of this charge density on the entire surface of the molecule will then generate a sigma profile (σ-profile) through the use of COSMOtherm software [52]. Then, the σ-profile will now be used as the basis by COSMO-RS to predict the desired thermodynamic properties. Nevertheless, it is to be noted that among the computational process being shown in Figure 2, a user is only required to insert the input, while all the computational process will be carried out by the software itself. Therefore, it is utmost important to input the right information to extract the desired output.

The input or simulation method of COSMO-RS in this work has been conducted by referring to the work of Kurnia et al. [39, 53]. Figure 3 shows the required input for calculating hydrogen bonding value before a proper simulation could be run.

Figure 3.

Inputs required to run an IL-hydrate system simulation in COSMO-RS.

As observed from Figure 3, the required inputs are temperature and the mole fraction of IL-hydrate system. For this work, the temperature is fixed at 10^°C, which is the normal temperature where hydrate will start to form. The effect of temperature is also proven not to be significant in this work, which will be explained later in the section of result and discussion. Next, the right value of mole fraction has to be entered for all four components including cation, anion, water, and involved gas. These mole fraction values need to be calculated beforehand as shown in Table 2. Similar to an experimental method that has been carried out by the chosen papers [25, 30, 49, 54], this simulation also considers that IL is inserted into the water at a mass fraction 10 wt%. Besides, since COSMO-RS considers IL is made up of equimolar cation and anion, a mole of IL will be divided equally into half a mole of cation and half a mole of the anion in the calculation [36, 55, 56].

Table 2.

Example of calculation of mole fraction.

^* Assuming 100 g of mixture and IL is inserted at a mass fraction of 10 wt%, then 90 g will be water and 10 g will be IL.

When all inputs are inserted, the simulation can now be run. Similar simulation method is applied for all other desired properties including sigma profile, activity coefficient, and solubility of IL in water. When all fundamental property value is collected, the next step is the identification of pattern and, later, the development of correlation using multiple regression analysis.

2.3 Prediction of inhibition ability of ammonium-based ILs

In total, 20 ammonium-based ILs have been selected for this study based on literature review. For cations, only shorter alkyl chains cations starting from tetramethylammonium up to tetrabutylammonium cations are chosen because longer cations are not effective [29, 30]. This might be because shorter alkyl chains are easier to be adsorbed by crystal surface. Longer alkyl chain, on the other hand, might even promote the formation of hydrates due to their increased hydrophobicity to react with water [57]. On the other hand, anions are made up of halide group and tetrafluoroborate [BF4]⁻ and hydroxide [OH]⁻ ions due to their strong electrostatic charges and tendency to form hydrogen bonding with water [30] (Table 3).

Table 3.

List of ammonium-based ILs being predicted.

All of the above chemicals will be simulated and calculated in COSMO-RS, which the calculations were carried out using TURBOMOLE6.1. The quantum chemical calculation follows the DFT, using the BP functional B88-86 with a TZVP basis set and the resolution of identity standard (RI) approximation.

3. Progress and discussion

3.1 Correlation development and validation

Using the four fundamental properties that have been identified earlier, an effort to relate them with the effectiveness of IL as a hydrate inhibitor has been carried out. These four properties are sigma profile, hydrogen bonding energy, activity coefficient, and solubility of IL in water. The following sections now thoroughly report and discuss if these four fundamental properties have successfully been related to IL inhibition ability.

3.1.1 Interpretation of sigma profile graphs

A sigma profile graph in COSMO-RS allows us to understand certain aspects of an IL-water system. The main information we can obtain from the graph is to learn about the hydrophobicity of IL and the tendency of IL to act as a hydrogen bond donor or hydrogen bond acceptor. According to Klamt [5, 58], the sigma profile graph can be divided into three regions. The first region is the hydrogen bond donor region (at the left of −1.0 e/nm²), the second region is nonpolar region (between −1.0 and 1.0 e/nm²), and the thirdly region is the acceptor region (at the right of 1.0 e/nm²). By judging at which region the peak of an IL locates, the tendency of IL to act as hydrogen bond donor or acceptor would be identified. Generally, a peak that locates at the right side of the sigma profile graph indicates the more electronegative area and acts as an H-bond acceptor.

Now, Figure 4 shows the sigma profile graph of EMIM-Cl, BMIM-Br, and water molecules. Looking at the sigma profile of water molecules as shown in Figure 4, it is observed that water has two high peaks, one in the hydrogen bond donor region and another in the acceptor region [39]. This indicates that water has a high affinity toward both acceptor and donor. Furthermore, Figure 4 shows the sigma profile of two ILs, which are EMIM-Cl and BMIM-Br. From the figure, it is observed that cations EMIM and BMIM both have their peak in the nonpolar region. However, water molecules which have peaks in the polar region tend to have higher affinity only with strong hydrogen bond donor or acceptor, but not cation that lays its peak in the nonpolar region [59]. As a result, cations do not interact much with water molecules. Meanwhile, anions that have their peaks in hydrogen bonding acceptor region are more attractive to water molecules. Hence, this inferred that anion is the main ion that interacts with water molecules to prevent hydrate formation, whereas cation merely contributes very slightly in the process [57].

Figure 4.

Sigma profile graph of EMIM-Cl and BMIM-Br.

Moreover, we can see that EMIM, which has a shorter alkyl chain length, has its peak nearer to the polar region than BMIM. As consequences, EMIM is also more polarized and hydrophilic than BMIM, which is a desired characteristic of a good hydrate inhibitor. This also proves that a cation with shorter alkyl chain length is preferable during the tuning of IL inhibitor, as a shorter cation is less bulky and hence can more effectively interact with water molecules [49, 60]. For anion, Cl⁻ proves itself to be a better H-bond acceptor as it has a peak at the right side of the graph, which is the indication of its further electronegative. This at the same time means that Cl⁻ will be more effective in accepting H-bond from water molecules than Br⁻. Therefore, this makes Cl⁻ more hydrophilic and serves as a better anion for hydrate inhibitor.

In short, referring to Figure 4, we can see that EMIM-Cl is the best combination of ions among the two types of ILs. Due to its lower alkyl chain length cation and a more electronegative anion, it should perform the best among the four ILs. This deduction is supported by the work of Xiao et al. [30], which reported the order of IL effectiveness as EMIM-Cl > EMIM-Br > BMIM-Cl > BMIM-Br.

3.1.2 Hydrogen bonding

Although hydrogen bonding strength has been widely quoted to have a relationship with the effectiveness of IL as hydrate inhibitor [29, 30], so far, no work has been conducted to prove this relationship. In this work, validation is done and has successfully proven that a linear relationship exists between hydrogen bonding strength and the effectiveness of IL as hydrate. This linearity is validated through four different sets of data that comes from three papers [25, 30, 61]. All four sets of data show good linearity relationship, with the highest regression value as R² = 1 and the lowest as R² = 0.8926. As a result, this implies that the prediction of IL effectiveness could be made through the comparison of hydrogen bonding strength.

Besides proving this relationship, several interesting findings have also been observed throughout the process. Firstly, computation of COSMO-RS, in total, will calculate three kinds of energy value for an IL, namely, misfit energy (E_MF), hydrogen bonding energy (E_HB), and van der Waals energy (E_vdW). The summation of these three energies leads to the value of total interaction energy (E_int). Although hydrogen bonding strength is known to affect the effectiveness of IL, the significance of other energies could not be neglected yet. Hence, in Figure 8, all types of predicted energies including E_MF, E_vdW, E_HB, and E_INT are plotted against average depression temperature to determine if these energies could also affect the effectiveness of ILs as hydrate inhibitor.

Figure 5 demonstrates that for ILs with BMIM cation, it is evidently shown the anion contributes more to the total interaction energy than the cation. The reason behind this is virtually consistent; van der Waals energies are nearly constant for all of the tested ILs and have thus no effect on the temperature depression. The contribution of misfit energy, having only a regression value of 0.2247, is also negligible. This leaves the hydrogen bonding energy to be the only energy that plays an essential role in affecting the effectiveness of BMIM-ILs. Furthermore, the relationship between total interaction energy (E_INT) and temperature depression is also not convincing. This graph hence supports the earlier statement that hydrogen bonding strength between cation and anion is the most important type of energy that regulates IL interaction with water molecules [29, 30, 62]. The same pattern of relationship is then also observed in another two data sets from the work of Xiao et al. [30]. Figure 6 now shows the relationship between average temperature depression of ILs and the predicted hydrogen bonding energy from COSMO-RS.

Figure 5.

Average temperature depression from Sabil et al. [25] work vs. types of predicted energy (binary components).

Figure 6.

Average temperature depression against predicted hydrogen bonding energy for both EMIM- and BMIM-based ILs (binary components).

Clearly, the graph shows that the temperature depression value of IL-hydrate system is directly proportional to the hydrogen bonding energy (E_HB) for both EMIM-based and BMIM-based ILs. The larger the absolute value of E_HB, the higher the temperature depression of a hydrate system. For instance, for BMIM-based ILs in this graph, the rank of E_HB from highest to lowest is as BMIM-Cl > BMIM-Br > BMIM-I > BMIM-BF₄. The same ranking occurred to the average temperature depression as well, where BMIM-Cl has the highest temperature depression and BMIM-BF₄ has the lowest depression. This ranking could be explained by the fact that among four anions, Cl⁻ anion has the highest polarized charge and thus acts as the best hydrogen bond acceptor. BF₄⁻ anion, on the other hand, has the lowest polarized charge after Br⁻ and I⁻ anion and thus shows the lowest hydrogen bond strength because it is the weakest hydrogen bond acceptor among all. This graph, however, also displays an interesting finding, which is the separation of EMIM- and BMIM-based ILs into two different data sets, instead of one. This step is necessary as the combination of all ILs into one data set may lower the linearity of relationship. This statement is supported by Figure 10, which shows a graph of average temperature depression against predicted hydrogen bonding energy.

Figure 7 inferred that linear relationship only exists when ILs with the same cation are compared. An early deduction is that to ensure a linear relationship for a set of data, only one single ion, which is either cation or anion, can vary, while another one must be fixed. The relationship could not be applied to predict ILs with different cations and anions. This deduction is supported by Figure 8, which shows the regression value between average depression temperature and hydrogen bonding strength for a set of ILs with different cations but same Cl⁻ anion.

Figure 7.

Average temperature depression against predicted hydrogen bonding energy for a single data set consisting of both EMIM- and BMIM-based ILs (binary components).

Figure 8.

Average temperature depression against predicted hydrogen bonding energy for ILs with Cl- as anion but different cations (binary components).

With the regression value as high as 0.8976 from Figure 8, this supports our deduction earlier, where one ion must be fixed and another one could be varied to see the relationship. Furthermore, it is noticeable that when ILs with fixed anion but different cations are measured, the relationship between hydrogen bonding strength and average depression temperature is inversely proportional as before. The higher the absolute value of hydrogen bonding energy, the lower the average temperature depression. This could be explained by the sigma potential graph that has been discussed earlier. Previously, it is explained that cations generally have their peaks located in the nonpolar region of sigma profile, which is from −1 to 1 e/nm².

On the other hand, water molecules show two high peaks, one at the region of hydrogen bond donor and another at hydrogen bond acceptor. As a result, water molecules tend to have higher affinity only with strong hydrogen bond donor or acceptor, but not cation that lays its peak in the nonpolar region. Hence, this inferred that anion is the main ion that interacts with water molecules to prevent hydrate formation, whereas cation merely contributes very slightly in the process [57]. Since cations have a low affinity with water molecules, this also indicates that most of the cations in water will continue to bond with anions. In that case, the excess hydrogen bonding energy provided by stronger cation (that has higher E_HB) is unnecessary. This stronger hydrogen bonding energy will be used by cation to bond with anion, thus reduces the number of anions that are free to interact with water molecules. As a consequence, it will bring about an inverse effect on average temperature depression and reduce the effectiveness of ILs as a hydrate inhibitor.

In short, linear relationship does exist between hydrogen bonding strength and the thermodynamic hydrate inhibition ability of an IL. For a set of ILs with fixed cation and different anions, stronger hydrogen bonding between ILs lead to higher average depression temperature of an IL-hydrate system. Vice versa, for a set of ILs with fixed anion but different cations, stronger E_HB produces lower depression temperature. A lower depression of temperature subsequently signifies that the IL is less capable of shifting the equilibrium curve and is thus a weaker THI inhibitor. Predicted hydrogen bonding energy computed by COSMO-RS through a binary system consisting of only cation and anion has thus proven to be useful in predicting the effectiveness of ILs as inhibitors.

However, the above method of computation in COSMO-RS involves only the interaction between cation and anion, and it does not represent the hydrate system fully. Thus, the second computation of the quaternary system containing cation, anion, water, and methane gas has been conducted. Similar graphs have been plotted to find out how consistent E_HB is in predicting the effectiveness of IL. Figure 9 shows the graph of average depression temperature plotted against a different type of predicted energies.

Figure 9.

Average temperature depression from Sabil et al. [25] work vs. types of predicted energy (quaternary components).

As observed from Figure 12, when quaternary components are involved, which include cations, anions, water, and methane, it is still obvious that hydrogen bonding energy (E_HB) is the main energy that influences the hydrate inhibition effect. Meanwhile, misfit energy and van der Waals energy have only a low regression value that is below 0.10. However, it is noticed that total interaction energy (E_INT) provides a slightly higher regression value than E_HB which is 0.6848 than 0.6671, which does not occur in binary component simulation. This could be because while involving more components such as methane and water, the van der Waals energy and misfit energy between different components are now more significant and influential. As compared to binary component regression value, the highest regression value that is obtained here is only 0.6848, which is extracted from the E_INT. Nevertheless, this low regression value could be improved to 0.8276 by removing the outlier which is BMIM-HSO₄ (1-butyl-3-methylimidazolium hydrogen sulfate) as shown in Figure 10. This is because of the nature of HSO₄⁻ anion, which has an extra hydrogen bonding functional group, OH⁻ (hydroxide), and thus resulting in stronger inhibition effect [62].

Figure 10.

Average temperature depression from Sabil et al. [25] work vs. predicted total interaction energy (quaternary components, without BMIM-HSO₄).

Figure 11 then shows the regression value of two more data sets from the work of Xiao et al. [30]. For both sets of data, total interaction energy (E_INT) gives the highest regression value too.

Figure 11.

Average temperature depression from Xiao et al. [30] work vs. predicted total interaction energy (quaternary components).

Similarly, from Figure 12, when the anions are fixed and cations are varied to study, the temperature depression value also decreases as the hydrogen bonding energy becomes more negative (stronger) which leads the increase in total interaction energy. Therefore, generally, COSMO-RS simulation of binary components and quaternary components both work well as a quick prediction for the effectiveness of IL as a thermodynamic hydrate inhibitor. Table 4 shows the regression value of both binary and quaternary components simulation.

Figure 12.

Average temperature depression against predicted total interaction energy for ILs with Cl- as anion but different cations (quaternary components).

Table 4.

Comparison of regression values produced by binary and quaternary components simulation.

Here, it is shown that the simulation of hydrogen bonding energy (E_HB) of binary components simulation provides a more consistent regression value. On the other hand, total interaction energy (E_INT) of quaternary components simulation more accurately reflects out the hydrated state which involves not only the IL itself but also water molecules and methane gas. To determine whether binary or quaternary components simulation is more effective in predicting ILs effectiveness, more sets of experimental data should be validated using the above approach. However, experimental work that tested ILs set with fixed anion or cation is very limited. Therefore, it is hard to conclude here whether binary or quaternary components simulation is more superior. Nevertheless, since real hydrate system consists of the interaction between water, methane, and IL, quaternary components simulation will be further studied, and correlation will be developed in this work.

From the previous analysis for quaternary simulation, it is observed that the total interaction energy of anion and cation has a different effect on average temperature depression. Anion with higher interaction energy shows a higher average temperature depression, while the stronger interaction energy of cation reduces the average temperature depression. Due to the opposite effect of these two types of interaction energies (cation and anion), it is thus a must to consider them separately during the development of correlation. This results in the splitting of total interaction energy (E_INT) into two variables, which are E_INT contributed by anion (E_{INT, A}) and E_INT contributed by cation (E_{INT, C}). Both of them are available and obtainable from COSMO-RS simulation. Table 5 shows an example of E_{INT, C}, E_{INT, A}, and E_INT calculated by COSMO-RS for the ILs from the work of Sabil et al. [25].

Table 5.

Type of interaction energies predicted by COSMO-RS for the work of Sabil et al. [25].

From Table 5, it is clear that the summation of E_{INT, A} and E_{INT, C} would result in the value of E_INT. In comparison, it is also evidently shown that anion contributes more to the total interaction energy than the cation. Now after obtaining the two variables, Minitab is used to assist in developing a suitable correlation for the prediction of average temperature reduced by each IL. Generally, the model could be described as.

$$Y={\mathrm{\beta }}_1+{\mathrm{\beta }}_2{X}_2+{\mathrm{\beta }}_3{X}_3+\dots$$ E2

$$\mathit{\Delta }\mathit{T}={\mathrm{\beta }}_1+{\mathrm{\beta }}_2{E}_{\mathrm{INT},\mathrm{C}}+{\mathrm{\beta }}_3{E}_{\mathrm{INT},\mathrm{A}}$$ E3

Among the many equations that have been tested, the best equation is listed below. It involves both E_{INT, A} and E_{INT, C} as independent variables.

Model:

$$\mathrm{\Delta }T=1.758+0.0643E_{\mathrm{INT},\mathrm{C}}-0.00559E_{\mathrm{INT},\mathrm{A}}$$ E4

Table 6 then shows the experimental value obtained from the literature review, as well as the predicted temperature using the above equation. It listed 25 ILs with their experimental average temperature depression value from 4 literature review [25, 30, 49, 54]. Using the values of ions’ interaction energy (E_INT,A and E_INT,C) obtained from COSMO-RS, average temperature depression has been predicted for each IL. Absolute error between the experimental and predicted value is then calculated and shown at the last row of the table, without considering the three extreme outliers that are highlighted in red.

Table 6.

Experimental and predicted average temperature depression of ILs for selected literature review.

Table 6 shows several interesting findings and limitations of the model. First, regarding the three extreme outliers, all three of them are substituted cations that have a hydroxyl (OH⁻) group. This type of substituted cation, as calculated by COSMO-RS, has an overly high E_{INT, C} (42.17 kJ/mol for [OH-C2MIM]-Cl as compared to 20.60 kJ/mol for EMIM-Cl), which is supposed to reduce their inhibition ability. But, in truth, hydroxyl group-substituted cation has constantly performed better than common cation because the OH⁻ serves as a strong hydrogen bond donor that will react with water [61, 63]. The increased interaction with water molecules will thus improve the average temperature depression [64]. Due to this reason, a large discrepancy is observed between experimental and predicted temperature depression for OH⁻-substituted cation-based ILs. This also signifies that the model developed earlier does not apply to hydroxyl group-substituted cations or possibly any other substituted cations ILs.

Next, a pure error which is caused by inconsistency between experiments has also limited the accuracy of this model. Table 7 shows the simplified list of ILs which have different experimental average temperature depression value obtained from the literature review.

Table 7.

ILs with inconsistent experimental average temperature depression.

As observed from Table 8, the experimental value obtained from literature review does not agree with each other. They are inconsistent, and this has thus hindered the development of a fully accurate model that could predict the inhibition ability of ILs as THI inhibitors. For instance, the inhibition ability of BMIM-BF4 was reported in three different papers, and the difference of experimental value from each paper is fairly large, ranging from 0.270 to 0.858°C. Nevertheless, Zare et al. [54] reported an experimental value of 0.460°C, which only presents a 2.92% error when compared to the predicted value.

Table 8.

Model summary and ANOVA for the model developed.

Next, looking at BMIM-HSO₄ and EMIM-HSO₄, it is experimentally proven that EMIM-HSO₄, which has a smaller alkyl chain length for cation, would serve as a better inhibitor [49, 60, 62]. However, because experimental values are obtained from two different papers [25, 54], BMIM-HSO₄ recorded a higher average temperature depression. This contradiction due to inconsistency again hardened the process of model development. Two factors could probably explain this inconsistency between experimental values: (i) purity of ILs being used in an experiment and (ii) experimental procedure and atmospheric condition.

In short, hydrogen bonding energy is the main type of energy that affects the interaction of ions with water and subsequently the inhibition ability of ILs. For a quaternary component simulation, however, total interaction energy shows a better linear relationship with average temperature depression. The model developed which considers cation interaction energy and anion interaction energy sufficiently predicts average temperature depression with an average error of 20.49%. It is to be noted that to a certain degree, the inconsistency between experimental values also contributed to the average error. Table 9 shows the regression statistics and P-value from the ANOVA test for the equation developed. The confidence level for the model is set at 95%, and thus a P-value of 0.000 (<0.05) signifies a reliable model.

Table 9.

Percentage difference of predicted ∆T for EMIM-Cl due to temperature difference.

3.1.3 Effect of temperature on predicted inhibition ability

From the earlier section, it is mentioned that the simulation work in this study is fixed at a temperature of 10°C, which is a common temperature where hydrates start to form. In this section, the effect of temperature is further examined to investigate if the predicted inhibition ability of ILs changes dramatically with temperature. Figure 13 shows the graph of predicted average temperature depression against simulation temperature.

Figure 13.

Graph of predicted average temperature depression against simulation temperature.

Nevertheless, if the percentage of difference is calculated out, it will be noticed that the effect of temperature is very insignificant. For example, for EMIM-Cl, using simulation of 10°C as reference state, the percentage difference for each temperature is shown at the table.

As shown in the table, the range of predicted average temperature depression is between 1.079 and 1.151°C, where the difference is really small. Furthermore, it is found out that most experimental studies involve only temperature range of −3.15 to 16.85°C (270–290 K) [3, 29, 30, 49, 63]. This means that the highest percentage difference is just around 3.12% (for −5°C). Hence, it can be concluded that the effect of temperature is insignificant and would not affect the screening process of ILs using the correlation.

3.1.4 Activity coefficient

As discussed by Kurnia et al. [39], the lower the activity coefficient of a water-IL mixture, the higher the interaction between components in the mixture. Khan et al. also explain that for a water-IL mixture, activity coefficient below 1 signifies favorable interaction between water and ILs in the mixture [34]. When ILs interact well with water, supposedly, less water will be free to bond with each other to form hydrate. Theoretically, the activity coefficient could then reflect out the inhibition ability of IL. Therefore, validation effort was made through four sets of data [25, 30, 61] to find out if the relationship between activity coefficient and average temperature depression exists. Figure 14 shows the graph of average temperature depression against the natural logarithm of activity coefficient.

Figure 14.

Graph of average temperature depression against ln(Y_W) for three data sets.

As shown in Figure 14, the highest regression value is observed for BMIM-based ILs from the work of Xiao et al., which is a mere 0.6658. Meanwhile, another two sets of data record unacceptably low regression value of only 0.0045 and 0.2989. Hence, regrettably, these three sets of data could not exhibit any significant relationship between these two variables. Nevertheless, a general pattern of decreasing average temperature depression is observed when the natural logarithm of activity coefficient increases (activity coefficient increases). The incapability of the activity coefficient in reflecting the inhibition ability of IL since the calculations of activity coefficient in COSMO-RS considers only the input of temperature, but no input of pressure is allowed. Meanwhile, in reality, hydrate occurs at low temperature (around 10°C) but high pressure. This kind of special nature of hydrate formation has thus made it hard for COSMO-RS to accurately predict out the activity coefficient water for a system of low temperature yet high pressure.

3.1.5 Solubility

A more soluble IL in water signifies that the IL can easily dissolve itself and interact with water molecules. Supposedly, a good IL should have high solubility in water, to bond with other water molecules and reduce the possibility of free water molecules from forming hydrate. To test the validity of the statement, four sets of data [25, 30, 61] were studied to find out if the relationship between the solubility of IL in water and average temperature depression exists. Figure 15 shows the graph of average temperature depression against the solubility of IL in water.

Figure 15.

Graph of average temperature depression against solubility of IL in water.

Here, the regression values for all three data sets are very low as well, with the lowest regression value of 0.057. Hence, similarly to the activity coefficient, no significant relationship could be deduced from this variable.

3.2 Prediction of inhibition ability of ammonium-based ILs

From the validation part, it has been confirmed that sigma profile and total interaction energy of ILs can be correlated to the effectiveness of an IL as THI inhibitor. Hence, in this section, prediction work will be conducted on 20 ammonium-based ILs (refer to Table 3) to determine their ability as hydrate inhibitor, through the study of their sigma profile and total interaction energies.

3.2.1 Sigma profile

Although sigma profile could not directly compute a value to represent the effectiveness of an IL as a hydrate inhibitor, it does show the affinity of an IL toward the water. The higher the affinity of IL toward the water, the more hydrophilic it is, and the easier it could interact with water. This will then result in a more effective hydrate inhibitor. Hence, in this section, three sigma profile graphs will be used to determine the affinity of each ammonium-based ILs toward the water. The first figure, Figure 16, displays the sigma profile of the four types of cations involved here, which range from TMA to TBA.

Figure 16.

Sigma profile of ammonium-based cations.

From Section 3.1.1, it is discussed that the sigma profile graphs could be divided into three regions: hydrogen bond donor region (at the left of −1.0 e/nm²), nonpolar region (between −1.0 and 1.0 e/nm²), and hydrogen bond acceptor region (at the right of 1.0 e/nm²). Here, all tetraalkylammonium-based cations have their peaks within the nonpolar region and are thus deduced to have a low affinity with water. This is because water molecules have only peaks within the hydrogen bond donor and hydrogen acceptor region. Due to this property, they do not interact well with ions that have a peak in the nonpolar region. However, when compared among themselves, TMA cation which has its highest peak at around −0.9 e/nm² performs the best because its peak is nearest to the polar region and thus has the highest affinity toward the water. This is because TMA has the lowest alkyl chain length, thus is less bulky and can easily interact with water molecules [34]. This makes TMA the most suitable cation among the four to be tuned as a hydrate inhibitor.

Figure 17 shows the sigma profile of five types of different anions. From this graph, it is observed that all anions have their peaks located in the polar region at the right side, which is the hydrogen bond acceptor region [33]. This indicates that all of them are electronegative and has a lone pair ready to share with another hydrogen bond donor. Due to their readiness to interact with hydrogen bond donor, they have high affinity with water molecules and tend to bond well with water molecules. The highest tendency of interaction goes to OH⁻ ion, which has its peak at 3.6 e/nm². In general, an anion that lays its peak further at the right side of the sigma profile graph is effective in inhibiting as it has high affinity with water molecules. This is because the further the peak to the right, the larger the sigma value and, thus, the more electronegative an anion is. The high electronegativity then results in higher interaction energy and thus interacts better with the water molecules. Meanwhile, BF₄⁻ ion that has its peak close to the nonpolar region is not an effective inhibitor anion because of its low polarized charge.

Figure 17.

Sigma profile of anions.

Lastly, the third figure, Figure 18, has selectively displayed the sigma profile graph for four ILs, including TMA-OH, TEA-Cl, TPA-I, and TBA-BF₄. The idea of this graph is to showcase several possible combinations of ILs by tuning the cation and anion. Here, it is easily observed that all cations show their peak in the nonpolar region. TMA cation shows its peak closest to the polar region and is thus the most suitable cation, due to its higher affinity with water. This could be explained by its short alkyl chain length as compared to others, which makes it more hydrophilic. Meanwhile, all anions lay in the polar region on the right side. The most electronegative anion is OH⁻ ion that has its peak furthest at the right. Due to its highest electronegativity and hence high interaction with water, it serves as the best anion to be used for an inhibitor. Therefore, from the graph, it is identifiable that TMA-OH is the best combination of all. This is followed by TEA-Cl, TPA-Br, and, finally, TBA-BF₄. From this graph, it is inferred that to choose the right anion for the hydrate inhibitor, its peak should be located as far as possible at the right side of the graph. This indicates a highly electronegative anion that can bond well with water molecules. Meanwhile, it is reported that most of the ILs cations have their peaks located in the nonpolar region. This characteristic causes cations to behave as nonpolar molecules that are hydrophobic and does not interact well with water molecules [58]. Therefore, a cation with the lowest hydrophobicity should be chosen to be tuned as a hydrate inhibitor, so that it will not hinder interactions between IL and water molecules. This, in turn, signifies that the most recommendable cation should have its peak closest to the left polar region.

Figure 18.

Sigma profile of several ammonium-based ILs.

3.2.2 Total interaction energies

In Section 3.1.2, a correlation has been developed to describe the relationship between the average depression temperature of IL-hydrate system and the total interaction energies. It is found out that both cation and anion interaction energy have a different effect on IL inhibition ability. High anion interaction energy is preferable, while high cation interaction energy will reduce an IL inhibition ability.

$$\mathrm{∆}T=1.758+0.0643E_{\mathrm{INT},\mathrm{C}}-0.00559E_{\mathrm{INT},\mathrm{A}}$$ E5

Using the above correlation, the ability of ammonium-based ILs has been predicted through the calculation of average temperature depression. Table 10 shows the list of ammonium-based ILs together with their total interaction energies and predicted inhibition ability measured in terms of average temperature depression.

Table 10.

Predicted average temperature depression of AILs.

Table 10 shows a list of tetraalkylammonium-based ILs, which range from cation tetramethylammonium to tetrabutylammonium paired with five types of different anions that are hydroxide ion, a tetrafluoroborate ion, chloride ion, bromide ion, and iodide ion. From this table, it is observed that when the anion is fixed, an increase in cation interaction energy, which is caused by the increase in alkyl chain length, will reduce average temperature depression. This again agrees to the earlier statement which explained that the longer alkyl chain length of cation, the bulkier it is and thus harder for it to interact with water molecules [49, 60]. This, as a result, increases its hydrophobicity, reduces its ability to bond with water, and is thus a less effective thermodynamic hydrate inhibitor [62]. In fact, among the five TBA ionic liquids (ILs), three of them show negative temperature depression. This is because of the poor combination of the bulky cation (TBA) and weak electronegativity anion (Br⁻, I⁻, BF₄⁻), resulting in a super ineffective inhibitor. A negative temperature depression signifies that instead of serving as hydrate inhibitor, they have now become hydrate promoter that favors the formation of the hydrate phase.

In terms of the effect of anion, we can see that the higher the interaction energy of anion, the higher the average temperature depression is. Here, the rank of E_INT is as OH⁻ > Cl⁻ > Br⁻ > l⁻ > BF₄⁻. This resulted in the average temperature depression to follow the same pattern. For instance, looking at tetramethylammonium ILs, inhibition ability rank is as TMA-OH > TMA-Cl > TMA-Br > TMA-I > TMA-BF₄. Hence, this again proves that the interaction energy provided by anion plays a crucial role in determining its inhibition ability. Also, this prediction agrees well with work reported by Tariq et al. [62]. In his review work, he reported that for a methylimidazolium-based IL, the order of efficiency follows as such C₂C₁im-Cl > C₂C₁im-Br > C₂C₁im-I >C₂C₁im-BF₄. Regrettably, OH⁻ ILs are not studied in Tariq’s work; yet, the whole ranking ranging from Cl⁻ to BF₄⁻ is similar to the predicted ranking. This proves that the developed correlation is performing outstandingly in predicting the inhibition ability of ILs. Lastly, from this model, TMA-OH is identified to show the strongest ability as THI, with the highest depression temperature of 1.97°C. This is due to the highly electronegative OH⁻ anion that bonds well with water molecules and a short alkyl chain length TMA cation that does not hinder the IL interaction with water molecules.

4. Conclusions

In conclusion, among the four identified fundamental properties, sigma profile and hydrogen bonding energy have been successfully correlated to the inhibition ability of IL. Sigma profile provides a qualitative understanding of each IL in the sense of their affinity toward water molecules. Meanwhile, hydrogen bonding energy, or later upgraded to total interaction energy, has been able to satisfactorily predict out a quantitative value of average temperature depression provided by each IL. This value will then tell us the effectiveness of each IL as a thermodynamic hydrate inhibitor. The correlation developed is validated with open literature and is found out to have an average error of 20.49%. From the predicted data, it is observed that TMA-OH depresses the temperature of IL-hydrate system by 1.97°C, whereas the widely studied EMIM-Cl can only experimentally depress the system by 1.22°C. TMA-OH has shown the highest inhibition ability due to the combination of its short alkyl chain length cation and a highly electronegative OH⁻ anion. Findings, however, show that this correlation is not suitable to be used for substituted cations, as the introduced functional group such as hydroxyl group will provide extra H bonding with water molecules. COSMO-RS simulation, on the other hand, has been proven to be applicable in computing fundamental properties of IL-hydrate system. Simulation of COSMO-RS in calculating fundamental properties paired with the correlation developed in this work could now serve as a pre-screening tool of ILs inhibition ability. This helps to narrow down the scope of ILs to be focused during experimental work and thus speeds up the rate of potential ILs being tested and applied to industrial processes.

Acknowledgments

YUTP financially supports this work under the Grant no. 015LCO-154.

Acronyms and abbreviations

ILs

ionic liquids

LDHI

low-dosage hydrate inhibitor

TMA-OH

tetramethylammonium hydroxide

TMA-Cl

tetramethylammonium chloride

TMA-Br

tetramethylammonium bromide

CDA

carbamyldicyanomethanide

N_2,2,2,2-Cl

tetraethylammonium chloride

BMIM-ClO₄

1-butyl-3-methylimidazolium perchlorate

BMIM-Cl

1-butyl-3-methylimidazolium chloride

BMIM-N(CN)₂

1-butyl-3-methylimidazolium dicyanamide

KHI

kinetic hydrate inhibitors

PVP

polyvinyl pyrrolidinium

∆T (°C)

average suppression temperature

COSMO-RS

conductor-like screening model for real solvent

CH₄

methane

BMIM-BF₄

1-butyl-3-methylimidazolium tetrafluoroborate

N_1,1,1,1-Cl

tetramethylammonium chloride

BMIM-CF₃SO₃

1-butyl-3-methylimidazolium trifluoromethanesulfonate

BMIM-CH₃SO₄

1-butyl-3-methylimidazolium methylsulfate