Molecular Interactions between Asphaltene and Surfactants in a Hydrocarbon Solvent: Application to Asphaltene Dispersion

Heavy oil and bitumen supply the vast majority of energy resources in Canada. Different methods can be implemented to produce oil from such unconventional resources. Surfactants are employed as an additive to water/steam to improve an injected fluid’s effectiveness and enhance oil recovery. One of the main fractions in bitumen is asphaltene, which is a non-symmetrical molecule. Studies of interactions between surfactants, anionic, and non-anionic, and asphaltene have been very limited in the literature. In this paper, we employed molecular dynamics (MD) simulation to theoretically focus on the interactions between surfactant molecules and different types of asphaltene molecules observed in real oil sands. Both non-anionic and anionic surfactants showed promising results in terms of dispersant efficiency; however, their performance depends on the asphaltene architecture. Moreover, a hydrogen/carbon (H/C) ratio of asphaltenes plays an inevitable role in asphaltene aggregation behavior. A higher H/C ratio resulted in decreasing asphaltene aggregation tendency. The results of these studies will give a deep understanding of the interactions between asphaltene and surfactant molecules.


Introduction
Precipitation of asphaltene has always been a detrimental process in the oil industry, and a massive amount of effort has been made to keep it dissolved in the oil phase [1][2][3]. Precipitation and aggregation of asphaltenes depend on various parameters, e.g., the thermodynamic condition of a system and architecture and structure of asphaltene molecules (non-symmetrical molecules) [4][5][6][7][8][9]. Asphaltenes are non-symmetrical molecules with wide ranges of sizes and structures, which make them hard to model. Various methods can be applied to prevent asphaltene precipitation during oil production, e.g., using inhibitors or dispersants [10]. Anionic and non-anionic surfactants can be used as inhibitors or dispersants [11][12][13][14][15][16]. In this regard, a thorough knowledge of the mechanisms and interactions between asphaltenes and surfactants plays a crucial role in formulating dispersants and inhibitors due to the complicated and broad ranges of asphaltene architectures. Such information also leads to a better selection of surfactant as a chemical-based enhanced oil recovery (EOR) [17][18][19][20][21].
Thanks to technology development, especially breakthroughs in computational capacity, molecular dynamics (MD) simulations draw tremendous attention. Via MD simulations, researchers have focused more precisely on different phenomena at the molecular scale, e.g., liquid-liquid and liquid-solid interactions [31][32][33] and interactions between fractions of heavy oils [34][35][36][37]. Several theoretical investigations have been performed to spotlight an aggregation process and colloidal behavior of asphaltene molecules in solvents [38][39][40][41][42][43][44][45][46]. Furthermore, asphaltene's structural effect on aggregation behavior has been studied in [47][48][49][50], and the molecular behavior of asphaltenes in their aggregates has also been investigated [35,51]. Tarefdar and Arisa [52] performed MD simulation to reveal thermodynamic effect on interactions between the resin and asphaltene. There are useful works in the open literature to find the mechanisms behind asphaltene aggregation [38,42,[53][54][55]. Headen et al. [56] used MD to examine the influence of resins on the colloidal behavior of asphaltene in different solvents, including toluene, n-heptane, and their mixtures. Tirjoo et al. [56] studied the effect of single-and divalent-ions on the aggregation behavior of asphaltene molecules in a solvent. They found that the presence of these ions is favorable for an asphaltene aggregation process. Celia-Silva et al. [57] used MD to address the aggregation behavior of asphaltene molecules in the presence of natural polymers as asphaltene dispersants. The dispersion behavior of these dispersants was favorable for both n-heptane and n-heptane-toluene mixtures; however, in toluene, the presence of these chemicals facilitated the formation of asphaltene dimers. In all cases, hydrogen bonding was the primary mechanism of the asphaltene dispersion process in the presence of chemicals.
Mizuhara et al. [58] investigated the effect of heteroatoms on asphaltene molecules' adsorption behavior at an interface of water and oil. Their theoretical outputs revealed that heteroatoms, including nitrogen, oxygen, and sulfur, facilitate the forming of hydrogen bonds between water molecules and asphaltenes. In addition, the highest van der Waals adsorption energy belongs to the structure with more sulfur content. Silva et al. [59] made an effort to address a fundamental question in breaking a water-oil emulsion by understanding how asphaltene, water, and emulsifier molecules interact with each other. Answering this question provides a better foundation and optimizes the required concentration of demulsifiers in a demulsification process. They concluded that an electric field does not significantly impact the demulsification process because using a low intensity electric field did not have dissolve salts in their simulation.
Moncayo-Riascos et al. [60] carried out both MD simulations and experiments to assess the influence of adding a methyl group to a solvent on rheological properties of asphaltenes. Based on their experimental and theoretical results, adding a methyl group resulted in significantly increasing the asphaltene-solvent solution's viscosity. The methyl group's position and quantity changed the asphaltene aggregates volume and shape, and these changes were the main contributors to this considerable difference in viscosity. Cao et al. [61] evaluated the Hildebrand solubility of maltenes and asphaltenes in pentane using MD simulations. They used different structures with different architectures to examine the structural effect of asphaltene on solubilization in pentane. Based on their simulation outputs, a higher amount of a benzene ring resulted in decreasing asphaltene solubility.
Furthermore, a larger aliphatic chain yields higher solubility in paraffinic solvents. Such outcomes could provide a better understanding for formulating a solvent for heavy oil or bitumen recovery, e.g., in the vapor extraction (VAPEX) method. Ahmadi and Chen [62] studied the interfacial behavior of surfactants, including all four types of surfactants, at an interface of water and asphaltenes. They used two asphaltene architectures, which were detected in Athabasca bitumen. Based on their outputs, anionic and non-anionic surfactants can reduce the probability of asphaltene aggregation; however, cationic and amphoteric surfactants cannot significantly affect the asphaltene aggregation.
No work has been done in the literature to focus on the effect of anionic and non-anionic surfactants on the asphaltene aggregation in the presence of a solvent, i.e., n-heptane, via molecular dynamics (MD) simulations. The main aim of the current work is to address this gap by employing MD simulations. This paper will provide a broader understanding of asphaltene aggregation in n-heptane in the presence of surfactants. Five different asphaltenes from different heavy oil samples [45,48,63] and two different surfactants, anionic (SDS) and non-anionic (TX-100), are used to address the aggregation trend of asphaltene in n-heptane very carefully. This paper's outcomes pave a road for formulating additives for using as a dispersant or inhibitor for asphaltene precipitation.

Force Field and Simulation Initialization
We used the Materials Studio software package (2020) [64] to carry out the molecular dynamics simulation. Table 1 provides the details, including parameters, ensembles, and methods, of the MD simulation. As reported in Table 1, the COMPASS force field was utilized in this study. The non-bonded interactions in this force field were Coloumbic and van der Waals interactions. The bonded interactions comprised different interactions, e.g., bond, angle, and torsion [65]. More information regarding the force field components and their definitions can be found in [37]. Figure 1 demonstrates the chemical formulas and structures of asphaltene molecules used in this paper. Table 2 reports details of all scenarios, including the simulation box's final size, system ID, surfactant, and asphaltene types and their H/C ratio. Figure 2 demonstrates the workflow of constructing a simulation box for performing MD simulations. As shown in Figure 2, five asphaltenes, five surfactants, and two hundred heptane molecules were placed in the simulation box with 40 Å × 40 Å × 40 Å dimension. As mentioned earlier, five different asphaltene molecules with different sources were used in the current work to evaluate the performance and efficacy of both anionic and non-anionic surfactants in the inhibition or dispersion of asphaltene aggregation in a heptane solution.    Figure 3 demonstrates the snapshot of the molecules' configuration in the simulation box for all fifteen systems at the equilibrium condition. Red molecules denote surfactants, blue molecules represent asphaltenes, and grey molecules stand for heptane. As demonstrated in Figure 1, the constructed simulation box needed to be optimized first. We performed the steepest descent algorithm for 200,000 iterations with convergence tolerance of displacement of 10 −1 Å and energy of 2 × 10 −5 kcal/mol to optimize the simulation box. After the geometry optimization process, we applied the NPT ensemble at a pressure of 1 Mpa and 298 K for 200 ps on each simulation box to equilibrate the simulation box. The configurations shown in Figure 3 were captured at the end of the equilibration process.    Figure 3 demonstrates the snapshot of the molecules' configuration in the simulation box for all fifteen systems at the equilibrium condition. Red molecules denote surfactants, blue molecules represent asphaltenes, and grey molecules stand for heptane. As demonstrated in Figure 1, the constructed simulation box needed to be optimized first. We performed the steepest descent algorithm for 200,000 iterations with convergence tolerance of displacement of 10 −1 Å and energy of 2 × 10 −5 kcal/mol to optimize the simulation box. After the geometry optimization process, we applied the NPT ensemble at a pressure of 1 Mpa and 298 K for 200 ps on each simulation box to equilibrate the simulation box. The configurations shown in Figure 3 were captured at the end of the equilibration process.  (e) (f)

Radial Distribution Function (RDF) Analysis
The RDF is defined as the ratio of the local density of molecules ρ(r) at distance r from a specific molecule at the origin to the average bulk density as follows [45]:

Radial Distribution Function (RDF) Analysis
The RDF is defined as the ratio of the local density of molecules ρ(r) at distance r from a specific molecule at the origin to the average bulk density as follows [45]: RDF analysis helps to have a better picture of the structural behavior of the specific molecule(s) inside a system. RDF plots help to reveal macro-molecular structures around specific molecules. A visible spike in an RDF plot means that a macro-molecular structure exists at that specific distance from the given atom or molecule. In addition, the area under the curve of an RDF plot proportionates to the cluster size around the given molecule. It means that the higher area under an RDF curve, the larger the cluster size. Furthermore, via RDF analysis along with visualization, we can understand the possible mechanisms which occur during the simulation time. Figure 4 depicts the RDF of asphaltene pair molecules for all systems. Based on the asphaltene molecular structure, different types of molecular interactions can occur. As illustrated in Figure 4a,b, the SDS surfactant had a lower RDF spike than the systems containing TX-100 and asphaltene without surfactant. However, in the cases of AS3, AS5, and AS6, no significant and meaningful difference between RDF plots for asphaltene in the presence of SDS and TX-100 could be observed (see Figure 4c-e). The probable reason behind this observation is because of the large asphaltene molecule size of these asphaltene molecules compared to AS1 and AS2. Hence, a macromolecular structure between asphaltene pairs in a short-range distance could not be observed for AS3, AS5, and AS6. Another reason to observe such results is because of a high H/C ratio (greater than one) of these asphaltene molecules, which means lower aromaticity. As reported in Table 2, only asphaltenes AS1 and AS2 had a H/C ratio of lower than one, revealing their higher aromaticity. The high aromaticity of asphaltene molecules favors for an asphaltene aggregation process. In other words, asphaltenes with high aromaticity are most likely to aggregate under specific thermodynamic conditions; however, asphaltenes with very low aromaticity probably do not tend to aggregate in wise ranges of thermodynamic conditions envelopes. That is why significant peaks could only have been observed for asphaltenes AS1 and AS2.
RDF analysis helps to have a better picture of the structural behavior of the specific molecule(s) inside a system. RDF plots help to reveal macro-molecular structures around specific molecules. A visible spike in an RDF plot means that a macro-molecular structure exists at that specific distance from the given atom or molecule. In addition, the area under the curve of an RDF plot proportionates to the cluster size around the given molecule. It means that the higher area under an RDF curve, the larger the cluster size. Furthermore, via RDF analysis along with visualization, we can understand the possible mechanisms which occur during the simulation time. Figure 4 depicts the RDF of asphaltene pair molecules for all systems. Based on the asphaltene molecular structure, different types of molecular interactions can occur. As illustrated in Figure 4a,b, the SDS surfactant had a lower RDF spike than the systems containing TX-100 and asphaltene without surfactant. However, in the cases of AS3, AS5, and AS6, no significant and meaningful difference between RDF plots for asphaltene in the presence of SDS and TX-100 could be observed (see Figure 4c-e). The probable reason behind this observation is because of the large asphaltene molecule size of these asphaltene molecules compared to AS1 and AS2. Hence, a macromolecular structure between asphaltene pairs in a short-range distance could not be observed for AS3, AS5, and AS6. Another reason to observe such results is because of a high H/C ratio (greater than one) of these asphaltene molecules, which means lower aromaticity. As reported in Table 2, only asphaltenes AS1 and AS2 had a H/C ratio of lower than one, revealing their higher aromaticity. The high aromaticity of asphaltene molecules favors for an asphaltene aggregation process. In other words, asphaltenes with high aromaticity are most likely to aggregate under specific thermodynamic conditions; however, asphaltenes with very low aromaticity probably do not tend to aggregate in wise ranges of thermodynamic conditions envelopes. That is why significant peaks could only have been observed for asphaltenes AS1 and AS2.    Figure 5 demonstrates the RDF of asphaltene-surfactant pair molecules for under-studied systems. As shown in Figure 5a, there was a spike around 1.5 Å in the case of SDS surfactant, which referred to the hydrogen bond between the SDS head group and the hydroxyl group of asphaltene A1. However, this did not occur for the TX-100 surfactant. Except for asphaltene AS2, all the rest of asphaltene molecules, including AS1, AS3, AS5, and AS6, had higher interactions with TX-100 due to the presence of a benzene ring in its structure.  Figure 5 demonstrates the RDF of asphaltene-surfactant pair molecules for under-studied systems. As shown in Figure 5a, there was a spike around 1.5 Å in the case of SDS surfactant, which referred to the hydrogen bond between the SDS head group and the hydroxyl group of asphaltene A1. However, this did not occur for the TX-100 surfactant. Except for asphaltene AS2, all the rest of asphaltene molecules, including AS1, AS3, AS5, and AS6, had higher interactions with TX-100 due to the presence of a benzene ring in its structure.   Figure 5 demonstrates the RDF of asphaltene-surfactant pair molecules for under-studied systems. As shown in Figure 5a, there was a spike around 1.5 Å in the case of SDS surfactant, which referred to the hydrogen bond between the SDS head group and the hydroxyl group of asphaltene A1. However, this did not occur for the TX-100 surfactant. Except for asphaltene AS2, all the rest of asphaltene molecules, including AS1, AS3, AS5, and AS6, had higher interactions with TX-100 due to the presence of a benzene ring in its structure.  Figure 6 demonstrates the final snapshots of systems at the end of the MD simulation; all heptane molecules are hidden for clarity purposes. As shown in Figure 6a, asphaltene AS1 molecules were away from each other and mostly interacted with SDS molecules, revealing SDS's good performance to disperse asphaltene molecules. In addition, as clearly seen from this figure, the SDS surfactant made a hydrogen bond with the hydroxyl group of asphaltene AS1; this observation is supported by the RDF plot in Figure 5a. However, as illustrated in Figure 6b, several asphaltene molecules tended to form an aggregate, which showed a lower efficiency of TX-100 in asphaltene dispersion. Figure  6c,d depict the final snapshots of molecules' configurations of systems AS2-SDS and AS2-TX, respectively. As illustrated in Figure 6c, most asphaltene molecules interacted with SDS and did not form any macro-molecular structure.
On the other hand, AS2 molecules formed a nanoaggregate, which shows weaker efficacy of TX-100 compared to SDS for dispersing asphaltene molecules. The story for AS3, AS5, and AS6 was different because there was no visible difference between snapshots captured for TX-100 and SDS for similar asphaltene. In other words, there was no tangible difference between TX-100 and SDS in asphaltene dispersion in cases having AS3, AS5, and AS6. However, partial interactions between SDS head and heteroatoms of asphaltene molecules for systems AS3-SDS, AS5-SDS, and AS6-SDS could be observed (see Figure 6e Figure 6 demonstrates the final snapshots of systems at the end of the MD simulation; all heptane molecules are hidden for clarity purposes. As shown in Figure 6a, asphaltene AS1 molecules were away from each other and mostly interacted with SDS molecules, revealing SDS's good performance to disperse asphaltene molecules. In addition, as clearly seen from this figure, the SDS surfactant made a hydrogen bond with the hydroxyl group of asphaltene AS1; this observation is supported by the RDF plot in Figure 5a. However, as illustrated in Figure 6b, several asphaltene molecules tended to form an aggregate, which showed a lower efficiency of TX-100 in asphaltene dispersion. Figure 6c,d depict the final snapshots of molecules' configurations of systems AS2-SDS and AS2-TX, respectively. As illustrated in Figure 6c, most asphaltene molecules interacted with SDS and did not form any macro-molecular structure.  Figure 6 demonstrates the final snapshots of systems at the end of the MD simulation; all heptane molecules are hidden for clarity purposes. As shown in Figure 6a, asphaltene AS1 molecules were away from each other and mostly interacted with SDS molecules, revealing SDS's good performance to disperse asphaltene molecules. In addition, as clearly seen from this figure, the SDS surfactant made a hydrogen bond with the hydroxyl group of asphaltene AS1; this observation is supported by the RDF plot in Figure 5a. However, as illustrated in Figure 6b, several asphaltene molecules tended to form an aggregate, which showed a lower efficiency of TX-100 in asphaltene dispersion. Figure  6c,d depict the final snapshots of molecules' configurations of systems AS2-SDS and AS2-TX, respectively. As illustrated in Figure 6c, most asphaltene molecules interacted with SDS and did not form any macro-molecular structure. On the other hand, AS2 molecules formed a nanoaggregate, which shows weaker efficacy of TX-100 compared to SDS for dispersing asphaltene molecules. The story for AS3, AS5, and AS6 was different because there was no visible difference between snapshots captured for TX-100 and SDS for similar asphaltene. In other words, there was no tangible difference between TX-100 and SDS in asphaltene dispersion in cases having AS3, AS5, and AS6. However, partial interactions between SDS head and heteroatoms of asphaltene molecules for systems AS3-SDS, AS5-SDS, and AS6-SDS could be observed (see Figure 6e Figure 7 shows the RDF of asphaltene-heptane pair molecules for all simulation boxes. As clearly seen from this figure, no significant spike could be observed, which means no significant cluster or aggregation occurs. However, in the case of having asphaltene AS1, AS2, and AS5, asphaltene molecules had greater interactions with n-heptane in the presence of SDS surfactant compared to TX-100 and no-surfactant cases, as depicted in Figure 7a-d. On the other hand, AS2 molecules formed a nanoaggregate, which shows weaker efficacy of TX-100 compared to SDS for dispersing asphaltene molecules. The story for AS3, AS5, and AS6 was different because there was no visible difference between snapshots captured for TX-100 and SDS for similar asphaltene. In other words, there was no tangible difference between TX-100 and SDS in asphaltene dispersion in cases having AS3, AS5, and AS6. However, partial interactions between SDS head and heteroatoms of asphaltene molecules for systems AS3-SDS, AS5-SDS, and AS6-SDS could be observed (see Figure 6e-i). Figure 7 shows the RDF of asphaltene-heptane pair molecules for all simulation boxes. As clearly seen from this figure, no significant spike could be observed, which means no significant cluster or aggregation occurs. However, in the case of having asphaltene AS1, AS2, and AS5, asphaltene molecules had greater interactions with n-heptane in the presence of SDS surfactant compared to TX-100 and no-surfactant cases, as depicted in Figure 7a-d.  Figure 7 shows the RDF of asphaltene-heptane pair molecules for all simulation boxes. As clearly seen from this figure, no significant spike could be observed, which means no significant cluster or aggregation occurs. However, in the case of having asphaltene AS1, AS2, and AS5, asphaltene molecules had greater interactions with n-heptane in the presence of SDS surfactant compared to TX-100 and no-surfactant cases, as depicted in Figure 7a-d.

Solubility Parameter
To have more comprehensive support for the results obtained by MD simulation, we performed a Monte Carlo simulation to obtain the solubility parameter of pure substances. The solubility parameter provides useful information regarding the state of the mixture in terms of solubility. In the present work, the solubility parameter helped us have a clear picture of the solubility of binary mixtures of asphaltene, surfactant, and heptane.
The cohesive energy density (CED) is defined as the required amount of energy for removing

Solubility Parameter
To have more comprehensive support for the results obtained by MD simulation, we performed a Monte Carlo simulation to obtain the solubility parameter of pure substances. The solubility parameter provides useful information regarding the state of the mixture in terms of solubility. In the present work, the solubility parameter helped us have a clear picture of the solubility of binary mixtures of asphaltene, surfactant, and heptane.
The cohesive energy density (CED) is defined as the required amount of energy for removing one molecule from solution to infinite separation. In 1963, Hildebrand proposed a parameter which is equal to the square root of the CED to quantify the solvency behavior of species [79]; the suggested parameter later became known as the "Hildebrand solubility parameter". The Hildebrand solubility parameter (δ) provides a numerical index for the solubility/miscibility state of components; species with similar solubility parameters are expected to be soluble. The solubility parameter (δ) of a specific molecule, e.g., asphaltene, surfactant, or heptane, can be determined by the following equation [80][81][82][83][84][85]: where ∆H v denotes the vaporization enthalpy, V m stands for the molar volume, R represents the gas constant, and T is temperature [85][86][87]. Table 3 reports the CED value and solubility parameters for asphaltenes, surfactants, and solvent at 1 MPa and 298 K. Forster et al. [88] proposed a strict threshold for solubility difference to see whether two components are miscible/soluble or not. According to this proposal, two substances are presumably miscible if their solubility difference (∆δ T ) < 2.0 MPa 1/2 . As reported in Table 3, for cases of AS1, AS2, and AS3, the difference between the asphaltenes' solubility parameters and the SDS solubility parameter was lower than the difference with a TX-100 solubility parameter. This shows a higher probability of having a miscible/soluble mixture in the case of having an SDS surfactant. However, the story was different for asphaltenes AS5 and AS6 due to a lower difference between their solubility parameters and the solubility parameter of TX-100.

Conclusions
The effect of both anionic and non-anionic surfactants on asphaltene precipitation in the heptane solvent was systematically evaluated using molecular dynamics simulations. In this regard, five asphaltene molecules with different architectures were used to include the potential impacts of an asphaltene structure on a precipitation process. According to the results achieved from the MD simulations, the following conclusions may be drawn:

•
The SDS surfactant had a lower RDF spike compared to the systems containing TX-100 and asphaltene without surfactant. However, for AS3, AS5, and AS6, no significant and meaningful difference between RDF plots for asphaltene in the presence of SDS and TX-100 could be observed. The probable reason behind this observation is because of the large asphaltene molecule size of these asphaltene molecules compared to AS1 and AS2. Hence, a macromolecular structure between asphaltene pairs in a short-range distance could not be observed for AS3, AS5, and AS6. Furthermore, asphaltene precipitation depends on asphaltene's aromaticity; the higher aromaticity, the higher tendency to be aggregated. In other words, a higher H/C ratio means lower aromaticity, which results in decreasing the aggregation tendency of asphaltene molecules. The aromaticity of asphaltenes AS1 and AS2 was higher than that of AS3, AS5, and AS6, which resulted in increasing the asphaltene aggregation behavior. • Both SDS and TX-100 surfactants could disperse both island (AS1) and archipelago (AS2) asphaltenes, revealing that they can be used as a dispersant for both types. Besides the architecture, the asphaltene molecular size and weight can play an essential role in forming a nano-aggregate with and without having surfactant dispersant.

•
Comparing solubility parameters of asphaltenes, surfactants, and heptane revealed that both surfactants were capable of being soluble into asphaltene aggregates. It means that these surfactants can disperse asphaltene molecules in the solvent and reduce their aggregation tendency. For asphaltenes AS1, AS2, and AS3, SDS surfactant was more likely to be miscible/soluble with asphaltenes; however, for asphaltenes AS5 and AS6, TX-100 was more soluble than SDS in an asphaltene-solvent solution.

Conflicts of Interest:
The authors declare no conflict of interest.