Energy Fuels 2009, 23, 5027–5035 Published on Web 09/11/2009
: DOI:10.1021/ef9004576
Aggregation and Partitioning of Model Asphaltenes at Toluene-Water Interfaces: Molecular Dynamics Simulations Tetyana Kuznicki,†,‡ Jacob H. Masliyah,† and Subir Bhattacharjee*,† †
Department of Chemical and Materials Engineering and ‡Department of Mechanical Engineering, University of Alberta, Edmonton, Alberta, T6G 2G6, Canada Received May 14, 2009. Revised Manuscript Received August 19, 2009
Molecular dynamics (MD) simulations were used to study the nanoaggregation of model asphaltene molecules in binary mixtures of toluene and water. Four types of model asphaltenes were studied: continental (C), Violanthrone-79 (VO-79), anionic continental (AC), and thiophenic anionic continental (TAC). Of these, the first two were uncharged, whereas the AC and TAC contained a side chain with a carboxylic terminal group. In all simulations, the model asphaltenes partitioned completely into the toluene phase of the phase-separated solvent mixture. The asphaltene molecules containing charged terminal groups remained tethered to the toluene-water interface (interface-bound), whereas uncharged asphaltenes aggregated within the central region of the toluene layer (core-bound). The stacked polyaromatic rings of the asphaltene aggregates showed a distinct preference to incline almost perpendicular to the toluene-water interface. Analysis of the mean square displacements of the asphaltenes reveals that for both charged and uncharged model asphaltenes, the lateral diffusion parallel to the toluene-water interface is dominant, whereas the transverse diffusion (normal to the interface) becomes significantly suppressed. This observation suggests that in thin films of toluene trapped between two aqueous phases, both interface-bound and core-bound asphaltenes, although subjected to vastly different interaction scenarios, exhibit similar diffusion behavior.
of these molecules or their aggregates with liquid-liquid interfaces, and their partitioning in liquid-liquid mixtures, is limited, and often empirical. This renders several petroleum processing operations to be based on essentially empirical foundations. For instance, oil-water emulsions encountered in petroleum processing are poorly understood owing to the presence of asphaltenes, with unpredictable, and often adverse, consequences on processes such as de-emulsification of oilwater emulsions, management of slop oil, choice of optimal diluents for bitumen and heavy oil, and disposal of wastes in deep wells. Considerable attention is being given toward understanding the molecular underpinnings of asphaltene dynamics and their aggregation. Since the recent review of Mullins,5 numerous analytical and experimental techniques were applied to study asphaltene nanoaggregation. Experimental techniques using nuclear magnetic resonance (NMR),7 centrifugation,8 direct current electrical conductivity,9 and small angle neutron scattering (SANS)10 revealed that the critical nanoaggregation concentration (CNAC) for different types of asphaltenes ranges between 50 and 150 mg/L and asphaltene behav;ior changes significantly at CNAC. Molecular dynamics simulations11,12 as
1. Introduction Asphaltenes are described as a solubility class of molecules that dissolve in excess toluene.1 Generally, they are a collection of diverse molecules consisting of polyaromatic units with varying amounts of aliphatic chains, few heteroatoms (such as, S, N, O), and traces of vanadium and nickel. Hydrogen/ carbon ratios of asphaltenes range between 1.0 and 1.2, with a heteroatom content of a few percent.2 They are ubiquitous in petroleum processing, where they are often cited as the “cholesterol” of petroleum,3 because they self-aggregate and eventually clog the production pipes and reservoir rocks. These compounds also play an important role in the global carbon cycle, as evidenced by the presence of derivatized, ionic asphaltenes in oceans in large quantities.4 Consequently, it is important to understand the physical and chemical properties of these entities to ensure efficient use and processing of hydrocarbon resources, as well as to assess the environmental fate of these molecules. Significant strides have been made toward understanding asphaltenes during the past decades.5 However, understanding of asphaltene aggregation in mixed solvents (for instance, aqueous-organic mixtures),6 interaction
(7) Lisitza, N. V.; Freed, D. E.; Sen, P. N.; Song, Y.-Q. Energy Fuels 2009, 23, 1189–1193. (8) Mostowfi, F.; Indo, K.; Mullins, O. C.; McFarlane, R. Energy Fuels 2009, 23, 1194–1200. (9) Zeng, H.; Song, Y.-Q.; Johnson, D. L.; Mullins, O. C. Energy Fuels 2009, 23, 1201–1208. (10) Headen, T. F.; Boek, E. S.; Stellbrink, J.; Scheven, U. M. Langmuir 2009, 25, 422–428. (11) Kuznicki, T.; Masliyah, J. H.; Bhattacharjee, S. Energy Fuels 2008, 22, 2379–2389. (12) Headen, T. F.; Boek, E. S.; Skipper, N. T. Energy Fuels 2009, 23, 1220–1229.
*To whom correspondence should be addressed. E-mail: Subir.B@ ualberta.ca. (1) Speight, J. G. The Chemistry and Technology of Petroleum [electronic resource], 4th ed.; CRC Press/Taylor & Francis: Boca Raton, 2007. (2) Vicente, L.; Soto, C.; Pacheco-Sanchez, H.; Hernandez-Trujillo, J.; Martinez-Magadan, J. M. Fluid Phase Equilib. 2006, 239, 100–106. (3) Creek, J. Energy Fuels 2005, 19, 1212–1224. (4) Dittmar, T.; Koch, B. P. Marine Chem. 2006, 102, 208–217. (5) Mullins, O. C. SPE J. 2008, 13, 48–57. (6) Khvostichenko, D. S.; Andersen, S. I. Energy Fuels 2008, 22, 3096– 3103. r 2009 American Chemical Society
5027
pubs.acs.org/EF
Energy Fuels 2009, 23, 5027–5035
: DOI:10.1021/ef9004576
Kuznicki et al. 13
well as the total combinatory Gibbs free energy techniques were applied to model asphaltenes in mixed solvents and observe asphaltene nanoaggregation there. These studies corroborate the promise of molecular modeling and simulations as powerful theoretical tools to address asphaltene behavior, as evidenced in some earlier studies.14,15 Large-scale molecular dynamics simulations of model asphaltenes in pure solvents (water, heptane, and toluene), as well as selected binary solvent systems, were recently reported.11 These types of simulations, however, are still considered to be at infancy in terms of capturing the entire breadth of asphaltene behavior. This is primarily owing to three critical factors inherent to simulations employing molecular models. First, there is no single molecule that can be considered as a suitable representation of asphaltenes.16-21 This results in a lack of consensus regarding whether single or multiple putative molecular structures can represent the collective behavior of asphaltenes. Second, the simulation box dimensions and time scales are often much smaller than the length and time scales over which meso-structures and macroscopic properties of asphaltenes and their aggregates evolve. Third, there are concerns regarding the robustness of the empirical force fields used in such simulations. Despite these limitations, the value of these simulations lies in their ability to provide some insight regarding the short-time dynamics, leading to improved understanding of the initial stages of the nanoaggregation process.7 The present study aims to test the applicability of these large-scale classical molecular dynamics (MD) simulations to describe the dynamics of nanoaggregation, and interfacial behaviors of the asphaltenes, and their aggregates at toluene-water interfaces. In this study, three putative molecular models of a continental asphaltene, with very similar molecular weights, similar hydrogen/carbon ratios, and identical heteroatoms (S, N, and O) were created. The basic continental model resembles one of the structures already presented in our previous study.11 Two variations of this basic continental structure were created, both of which had a terminal anionic carboxyl group on an aliphatic side chain. One of these ionic continental models was created with a thioether on a side chain, whereas another had a thiophenic ring. We also modeled the behavior of a commonly used surrogate of continental asphaltenes, namely, Violanthrone-79 (VO-79). Notably, VO-79, being a real molecule, can be used for experimental corroboration of the simulation results. These model molecules (continental, anionic continental, thiophenic anionic continental, and Violanthrone-79) were added to toluene-water mixtures, and their aggregation, diffusion, and interfacial behavior were studied employing MD simulations. A key feature of these simulations is the consideration of
large numbers of molecules, and use of fairly large systems, such that statistical information about the aggregates could be reliable. The simulation results were used to illustrate the combined influence of interfacial affinities and self-aggregation tendencies of the different model asphaltenes on the aggregate structures and the dynamics of these aggregates in the thin toluene film. 2. Molecular Dynamics Simulation Molecular dynamics (MD) simulations were conducted using the freeware GROMACS (version 3.3).22,23 For all MD simulations, the GROMOS96 force field24 was used. General aspects of adapting the software for simulating the petroleumlike systems were discussed in our earlier papers.11,25 Here, we present simulations of partitioning, aggregation, and interfacial dynamics of several model asphaltenes in toluene-water mixtures. A total of five simulations were conducted. Of these, three simulations were performed using putative molecules having the same continental type polyaromatic core, with variations in the surfactant-like and polar structural elements. These variations include anionic carboxylate terminal groups tethered to the aliphatic chains, and inclusion of sulfur as either thioether (in side chain) or thiophene (in polyaromatic ring). Two additional simulations were conducted using a molecular model of the well-known asphaltene surrogate, VO-79. These simulations explored the influence of VO-79 concentration on aggregation and interfacial dynamics. These two simulations also served as a benchmark, against which the behavior of the other putative asphaltene models were compared. 2.1. Simulation Procedure. The simulations were performed in the isothermal-isobaric (NPT) ensemble at temperature T = 298 K and pressure P = 105 Pa (1 bar). Weak coupling of the molecules to a solvent bath of constant temperature was maintained using the Berendsen thermostat with a coupling constant of τT = 0.1 ps.26 The pressure was controlled using the Berendsen algorithm with a pressure coupling constant of τp =1 ps, applying a compressibility of 4.510-5 kJ-1 mol nm3. Unless otherwise stated, the simulation parameters were generally maintained at the default values offered in GROMACS.27 A time step of 2 fs was used, and trajectories were recorded at 10 ps intervals in all simulations. The nonbonded potential truncation was performed using spherical cutoff method for the LJ interactions, with the cutoff separation fixed at 1.0 nm. The Particle Mesh Ewald (PME) technique was used (with a real space cutoff 1.0 nm) to handle the long-range electrostatic interactions. The Verlet neighbor list cutoff was 1.0 nm. It should be noted that the initial box dimensions are not important in these simulations, since we employ temperature and pressure coupling, which modify the box volume to attain the appropriate densities and compressibilities of the solvents. Trajectories were used to
(13) Johansson, B.; Friman, R.; Hakanpaa-Laitinen, H.; Rosenholm, J. B. Adv. Colloid Interface Sci. 2009, 147-48, 132–143. (14) Murgich, J.; Abanero, J.; Strausz, O. Energy Fuels 1999, 13, 278– 286. (15) Zhang, L. Q.; Greenfield, M. L. J. Chem. Phys. 2007, 127, 194502. (16) Murgich, J. Mol. Sim. 2003, 29, 451–461. (17) Strausz, O. P.; Safarik, I.; Lown, E. M.; Morales-Izquierdo, A. Energy Fuels 2008, 22, 1156–1166. (18) Mullins, O. C.; Martinez-Haya, B.; Marshall, A. G. Energy Fuels 2008, 22, 1765–1773. (19) Becker, C.; Qian, K.; Russell, D. H. Anal. Chem. 2008, 80, 8592– 8597. (20) Pomerantz, A. E.; Hammond, M. R.; Morrow, A. L.; Mullins, O. C.; Zare, R. N. Energy Fuels 2009, 23, 1162–1168. (21) Boek, E. S.; Yakovlev, D. S.; Headen, T. F. Energy Fuels 2009, 23, 1209–1219.
(22) Berendsen, H. J. C.; Vanderspoel, D.; Vandrunen, R. Comput. Phys. Commun. 1995, 91, 43–56. (23) Lindahl, E.; Hess, B.; van der Spoel, D. J. Mol. Model. 2001, 7, 306–317. (24) Gunsteren, V. The GROMOS home page: GROMOS96 Manual.; 2007. (25) Kuznicki, T.; Masliyah, J. H.; Bhattacharjee, S. Langmuir 2007, 23, 1792–1803. (26) Berendsen, H. J. C.; Postma, J. P. M.; Vangunsteren, W. F.; Dinola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684–3690. (27) Spoel, D. V.; vanBuuren, A. R.; Apol, E.; Meulenhoff, P. J.; Tieleman, D. P.; Sijbers, A.; Feenstra, K. The GROMACS Manual (V3.3); Groningen University: Groningen, The Netherlands, 2004.
5028
Energy Fuels 2009, 23, 5027–5035
: DOI:10.1021/ef9004576
Kuznicki et al.
compute the density distributions, system energy, radial distribution functions, and mean square displacements of different molecules. These were computed using appropriate postprocessing programs available in GROMACS,27 as well as using the visualization program VMD.28 2.2. Molecular Models and Simulation Box Construction. Four model molecules representing some key structural and compositional features of asphaltenes of molecular weight of approximately 800 were used in these simulations. These molecules are termed “continental” (C54H65NO2S, MW: 792), “anionic continental” (AC, C54H62NO4S, MW: 821), “thiophenic anionic continental” (TAC, C53H62NO4S, MW: 809) and “Violanthrone-79” (VO-79, C50H48O4, MW: 713). The two-dimensional structures of these molecules are shown in Figure 1. All four of the asphaltene models have one large aromatic region per molecule and several side chains. The anionic continental (AC) and the thiophenic anionic continental (TAC) molecules contain a terminal COO- group on one of the side chains. Furthermore, the S atoms are present in AC as thioether and in TAC as thiophene. VO-79 was used as one of the asphaltene model molecules that has some features of the C7 Athabasca asphaltenes.29 Three-dimensional structures of these model molecules were created and geometry optimized using Accelrys Materials Studio molecular modeling software. The coordinates were then exported and submitted to PRODRG30 to generate the topology files necessary for MD simulations in GROMACS. The alkane chain of the molecules was modeled as united atom structures (no explicit hydrogen atoms are modeled).31 All the aromatic regions were modeled by considering the sp2 hybridized carbons. For the anionic asphaltenes, one sodium ion per molecule was included in the molecular topology to render the simulation boxes electroneutral. Further details of the simulations of binary mixtures of water and toluene are available elsewhere.11 In all simulations, the total number of model asphaltene molecules was kept constant at 24. The concentration (molarity) of the asphaltene molecules, as well as the ratio of toluene and water, were kept almost constant in all the simulations. A summary of the simulations conducted, with pertinent model parameters is given in Table 1. The simulations were performed using water and toluene as the solvents. The simulation boxes were constructed by first packing an ordered array of 24 molecules of a single type of asphaltene molecules in a regular cubic lattice with an edge dimension of ≈8 nm. Following this, the box was filled with toluene. The box was then expanded in the x- direction to a length of ≈12 nm. The entire box was then filled randomly with water. This box construction is referred to as “dispersed in mixture” as the asphaltene molecules in this geometry are initially in random contact with both toluene and water molecules. We note that in all simulations, the number of toluene and water molecules was kept equal at 1405 and 12 946, respectively. During the model construction and subsequent simulations, 3D periodic boundary conditions were applied.
Figure 1. Models of asphaltene molecules used in this study. The C, AC, and TAC models (a-c) are putative. However, the VO-79 molecule (d) is a commercially available chemical, often employed as a synthetic surrogate of asphaltenes.
Details of simulation box construction for the binary solvent systems employed in this study are available elsewhere.25 The system energy and density were monitored during the simulations to ensure equilibration of the box. 3. Results and Discussion The primary objective of this study is to compare the aggregation dynamics of the four types of asphaltene molecules in the toluene-water system. The mixed solvent contains about 40% toluene by volume (mole fraction ca. 0.1), ensuring that the toluene phase (after solvent phase separation) will have sufficient volume to accommodate all the asphaltene molecules. The large simulation box size (and therefore, the large number of solvent molecules) was chosen to primarily accommodate this feature, and prevent any artificial box dimension induced partitioning of the asphaltenes. The asphaltene concentrations used in these simulations are above the reported CNAC of these compounds in this MW (700 800) range (simulations 1-4 of Table 1).7 3.1. Partitioning and Aggregation in Toluene-Water Systems. We first focus on the qualitative aspects of the asphaltene partitioning and aggregation behavior in toluenewater mixtures. Figure 2 depicts snapshots of the aggregate nanostructures formed in toluene-water mixtures at 7 ns of simulations 1-4 (Table 1). Since the snapshots are all at the same instant, visual comparison of the extent of asphaltene aggregation can be made between these four simulations for the different types of molecules. Although 7 ns is a very short time, and the snapshots merely depict the structure resulting from only the very high frequency (shortest time and length scale) dynamics, the images provide considerable insight regarding the aggregation and partitioning of these molecules in binary toluene-water systems. In Figure 2, the toluene molecules are intentionally not displayed so as to highlight the partitioning of the asphaltenes. Two features are common in
(28) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996, 14, 33. (29) Lopez-Linares, F.; Carbognani, L.; Gonzalez, M. F.; Sosa-Stull, C.; Figueras, M.; Pereira-Almao, P. Energy Fuels 2006, 20, 2748–2750. (30) Schuttelkopf, A. W.; van Aalten, D. M. F. Acta Crystallogr. Sect. D-Biol. Crystallogr. 2004, 60, 1355-1363, Part 8. (31) Sundararajan, V. Computational Modeling of Membrane Bilayers; Academic Press: London; Burlington, Mass., 2008.
5029
Energy Fuels 2009, 23, 5027–5035
: DOI:10.1021/ef9004576
Kuznicki et al.
Table 1. Different Types of Simulations Conducted in This Studya simulation
asphalt model
1 2 3 4 5
C AC TAC VO-79 VO-79
Nasphalt
time (ns)
final size (nm)
Dasphalt (m2/s)
24 24 24 24 5
10 9 7 20 20
11.9 7.6 7.6 12.2 7.5 7.5 12.2 7.5 7.5 12.1 7.5 7.5 11.5 7.6 7.6
(5.2 ( 1.0) 10-10 (4.8 ( 1.7) 10-10 (4.0 ( 0.6) 10-10 (4.7 ( 1.5) 10-10 (3.8 ( 2.0) 10-10
a The number of water and toluene molecules was constant in every simulation and equal to 12 946 and 1405, respectively. Diffusion coefficients of asphaltenes (Dasphalt) are averaged over the last nanosecond of the simulations.
aggregate structures of the C and VO-79 molecules might be attributed to their different molecular structures, notably, the polar nature of their polyaromatic ring regions (cf. Figure 1). It appears that the influence of the charge on the terminal carboxylic anions, located on the alkane side-chains of the model AC and TAC asphaltenes, overcomes the effect of the polarity of their polyaromatic fused planes. In other words, uncharged asphaltenes remain core-bound within the toluene film, whereas charged asphaltenes remain interfacebound. After the initial capture of carboxyl groups at the interface, stacking of polyaromatic rings of different model asphaltenes results in the formation of a complex 3D aggregate structure. It is possible that, in real systems, such an associated structure can easily provide considerable rigidity to the oil-water interface. In all these simulations (1-4), the initial relaxation of the total system energy (not shown) occurred within the first nanosecond of the simulation. The solvents demixed spontaneously within the first two nanoseconds to yield alternating layers of toluene and water. The partitioning of the asphaltenes into the toluene phase was also rapid, and was complete within about 2 ns. In other words, within 2 ns, the systems in simulations 1-4 all evolved from their starting configurations to provide the structural features closely resembling the snapshots of Figure 2. 3.2. Evolution of Larger Aggregates. Figure 3 depicts the VO-79 aggregate structures formed in toluene-water mixtures at 20 ns of simulations 4 and 5. In Figure 3a, the VO-79 asphaltene aggregate is suspended completely inside the toluene layer. Comparing this snapshot with the snapshot at 7 ns in Figure 2d, it is evident that between 7 and 20 ns, the entire population of VO-79 molecules self-associate to form a continuous rod-like aggregate. This growth of the aggregate into a single rod-like structure is representative of a long-term ordering behavior involving asphaltenes. It is, however, interesting to note that the VO-79 molecules did not associate with the toluene-water interfaces even at 20 ns. This lack if interfacial activity was also exhibited by the C asphaltenes in our previous work.11 Figure 3 b shows the aggregate structure of the VO-79 asphaltenes in simulation 5, where all conditions are identical to simulation 4, except that the concentration of VO-79 is approximately five times lower. Initially, all 5 asphaltene molecules were spread evenly over the entire volume of the initial simulation box. Over the 20 ns duration of the simulation, the asphaltene molecules sporadically aggregated, typically forming dimers and occasionally forming a trimer. No persistent aggregates are observed in this case, and majority of the asphaltenes kept migrating in isolation. All the asphaltene molecules remained within the toluene film during the simulation (after formation of the toluene layer). The aggregate morphologies indicate that asphaltene nanoaggregation is not observed below a CNAC. The system in simulation 5 was below the CNAC of VO-79 and did not
Figure 2. Aggregation and interfacial affinity of four model asphaltene molecules in toluene-water mixtures based on simulations 1-4 of Table 1. The snapshots in parts a-d correspond to the four different asphaltene models, namely, C, AC, TAC, and VO-79, respectively. All snapshots are taken at 7 ns. In all systems, the total number of molecules are identical. The toluene molecules are not shown to highlight the aggregate structure.
all simulation snapshots: (i) the toluene and water phase separate into two layers, and (ii) 100% of the asphaltenes partition into the toluene phase. Each type of model asphaltene showed slightly different aggregate nanostructuring. In Figure 2a, the aggregates of the continental molecules (simulation 1) remain suspended in the central region of the toluene layer, with some molecules still partly associated with the toluene-water interface. Out of 24 molecules, 10 molecules aggregated into two nanoaggregates containing 5 molecules each, and 6 molecules split evenly yielding 2 other nanostructures. Within the toluene layer, these molecules aggregated through stacking of their polyaromatic rings. These aggregates remained suspended in the central region of the toluene layer (core-bound) and did not associate with the toluene-water interface. The anionic continental molecules containing terminal carboxylic groups on one of the side chains (AC, Figure 1 b) exhibit a dramatically different partitioning and aggregation behavior. In this case, we observe in Figure 2b that nearly all of the AC molecules associate with the toluenewater interface, tethered through the deprotonated carboxyl groups. The same behavior is observed for the TAC asphaltenes (Figure 2 c). For both AC and TAC models, the anionic terminal groups act as tethers that bring the asphaltene aggregates to the toluene-water interfaces (interface-bound). However, the polyaromatic rings of these molecules still associate through π- stacking, and form clusters. Figure 2 d depicts the aggregation of VO-79. The behavior of VO-79 is quite analogous to that of the continental asphaltenes (Figure 2a), with the aggregates and their side chains remaining well-solubilized within the toluene layer (core-bound). The slight qualitative differences between the 5030
Energy Fuels 2009, 23, 5027–5035
: DOI:10.1021/ef9004576
Kuznicki et al.
exhibit any aggregation by sustained stacking of more than three polyaromatic rings. 3.3. Structural Analysis: Density Distribution and Pair Correlations. A more rigorous analysis of the time average density distribution and the aggregate nanostructures relative to the toluene-water interface in presented here. Density profiles obtained as a time average over the 6-7 ns duration of simulations 1-4 are shown in Figure 4. We considered the toluene-water interfaces in these density maps to be located where the density profiles of water (dashed line) and toluene (solid line) intersect. In Figure 4a and 4d, the C and VO-79 asphaltenes remain in the central region of the toluene layer, with the toluene molecules showing a greater density near the edges of the film. This confirms that uncharged asphaltenes are preferentially positioned in the center of the toluene film (core-bound), pushing the aromatic solvent to the interface. In contrast, both AC and TAC molecules concentrate near the interface (interfacebound) as shown in Figure 4b and 4c. It is interesting to note that the inclusion of sulfur (as thiophene or thioether) in the molecular structure does not influence the relative position of the model asphaltenes at the toluene-water interfaces, since there is no apparent difference in density profiles of AC and TAC molecules. Similarly, presence of NH does not result in significant deviation in the density profiles between the C asphaltene and VO-79 (which lacks the NH group in the polyaromatic ring). Figure 5 shows the pair correlation functions between several atomic groups of interest. The pair correlation function between the aromatic planes of the asphaltene models and toluene is shown in Figure 5a. Both continental and VO79 models of asphaltenes show very similar pair correlations with toluene at short-range. The slight differences in the
Figure 3. Concentration dependence of the aggregation behavior of VO-79 asphaltene molecules in toluene-water mixtures. (a) Snapshot of simulation 4 at 20 ns (containing 24 VO-79 molecules). (b) Snapshot of simulation 5 at 20 ns (containing 5 VO-79 molecules). The toluene molecules are not shown to highlight the aggregate structure. The number of toluene and water molecules in the two simulations are identical.
Figure 4. Density profiles of the different constituents of the simulated system plotted along the x-direction. The density profiles are averaged over 1 ns of the simulation time (between 6 and 7 ns): Parts (a) through (d) represent the results for simulations 1-4, respectively.
5031
Energy Fuels 2009, 23, 5027–5035
: DOI:10.1021/ef9004576
Kuznicki et al.
Figure 5. Pair correlation functions between various groups of atoms: (a) g(r) between the polyaromatic ring planes and toluene (simulations 1-4), depicting the slight structural differences induced by the continental and VO-79 molecules; (b) g(r) between the aromatic ring planes of C, AC, TAC, and VO-79 molecules and water; (c) g(r) between the deprotonated oxygen of terminal carboxyl groups and hydrogen of water (simulation 2); (d) g(r) between sulfur atoms of AC (thioether, simulation 2) and TAC (thiophene, simulation 3).
long-range radial distribution functions (RDF) may be attributed to the different polarities of the aromatic planes of VO-79 and C asphaltenes. The stronger polarity of the VO-79 aromatic region may induce longer range structures in the toluene molecules. The aromatic structures are the largest components of asphaltenes, and the interaction between aromatic planes occurs through π-π interactions. Because of this, both continental and VO-79 asphaltene models stack in toluene as observed in snapshots of simulations 1-4 (see Figures 2a and 3a). The stacked polyaromatic rings are typically at a separation of ca. 0.4 nm as noted in our previous work.11 It is interesting to note from Figure 5a that the toluene molecules are excluded from within the stacked polyaromatic platelets (the first peak of the radial distribution function appears at about 0.5 nm, which is greater than the stacking distance). Thus, although toluene is considered as the solubilizing agent for asphaltenes, our simulations indicate that these model asphaltenes aggregate by excluding toluene from between their polyaromatic rings. Figure 5b shows the pair correlation functions of the asphaltene polyaromatic planes with water for simulations 1-4. Continental and VO-79 molecules show a similar behavior, whereas the AC and TAC show similar trends. There are, however, marked differences between the uncharged (C and VO-79) and the charged (AC and TAC) asphaltenes with respect to their affinity for the toluenewater interface. The charged asphaltenes show a much higher affinity to the interface compared to C and VO-79 asphaltenes. Figure 5c depicts the pair correlation function between the deprotonated carboxyl group of the anionic C asphaltene and the hydrogen atoms of the water molecules. It is evident
that the deprotonated oxygen of the carboxyl groups of AC asphaltenes has a strong affinity toward the hydrogen of water. The sharp peak in g(r) at a distance of about 0.190.2 nm indicates that there is most likely a hydrogen bond forming between the deprotonated oxygen of COO- and hydrogen of water. The location of the peak corresponds to the literature reported value of the hydrogen bond length (approximately 0.197 nm).32 This hydrogen bonding competes with the π-π interactions of the aromatic planes and results in the interfacial activity of AC and TAC asphaltenes. Figure 5d depicts the sulfur-sulfur pair correlation function for the thioether asphaltene (AC, dashed line) and thiophenic asphaltene (TAC, solid line). The peaks in g(r) indicate that the sulfur atoms have a tendency to attract each other in both cases. The thiophenic sulfur atoms demonstrate higher probability to position close to each other compared to the thioether sulfur. This may be attributed to the fact that aliphatic chains with thioether have greater flexibility in toluene compared to the rigid polyaromatic ring structures containing thiophenes. From a comparison of these two RDFs, we observe that the long-range interactions of sulfur atoms are similar in AC and TAC. However, at short-range, thiophenes are more attracted to each other compared to thioethers. The RDF of thiophenes have more pronounced peaks and depressions within the first nanometer, indicating that the positioning of these atoms correspond strongly to the stacking distance of the polyaromatic rings in the aggregates. (32) Jeffrey, G. A. An Introduction to Hydrogen Bonding; Oxford University Press: New York, 1997.
5032
Energy Fuels 2009, 23, 5027–5035
: DOI:10.1021/ef9004576
Kuznicki et al.
the vector nB). Conversely, when sin(R) f 0, the aromatic planes are parallel to the interface. Figure 6b depicts the distribution of C, AC, TAC, and VO79 aromatic planes by plotting the distribution of sin(R) as a histogram. On the basis of the asphaltene distribution, we can conclude that all four types of asphaltenes (C, AC, TAC, and VO-79) prefer to orient themselves with high values of the angle of inclination relative to the toluene-water interface, that is, almost perpendicular to the interface. We denote this as the side-on configuration. A closer inspection of the histograms indicates that the continental model asphaltenes followed by VO-79 have the greatest tendency to organize perpendicular to the interface. This is attributed to the predominance of the π-stacking between the polyaromatic rings of these asphaltenes. The AC and TAC asphaltenes, on the other hand, exhibit a slightly broader distribution of angles, which might be an outcome of the competing influence of hydrogen bonding through the tethered anionic groups of these molecules or the sulfur atoms. The range of angles obtained in our simulations are in a good agreement with those obtained by Zhang and Greenfield33 in their simulations. We note that some molecules showed very slow rotation of their aromatic planes with respect to the toluene-water interface while some moved significantly over the duration of the simulations. Table 2 shows the mean values of sin(R) obtained from analysis of the histograms in Figure 6b. These were calculated assuming both normal and log-normal distributions. Both distributions yield identical mean values. For both the C and VO-79 molecules, the mean value of sin(R) was approximately 0.8. This is partially due to the fact that both molecules are uncharged, and exhibit the same tendency to stack inside the toluene layer, eventually yielding a large rodlike aggregate structure aligned parallel to the interface. This causes the individual polyaromatic rings to remain nearly perpendicular (side-on) to the interface. For AC and TAC molecules, we observe a lower value of sin(R), which is closer to 0.7. Although this still implies that the polyaromatic rings of these ionic molecules prefer to align perpendicular to the interface (side-on configuration), their preferential tethering to the oil-water interface prevents a large number of these molecules to stack efficiently. Consequently, a larger standard error or variance is observed for these molecules. 3.5. Mean Square Displacement and Diffusion. The mean square displacement trajectories of the molecules provide a facile understanding of the molecular dynamics of the various components of the multiphase system, and provides the self-diffusion coefficient of the molecular species. The selfdiffusivity of the molecules was calculated from the Einstein relation, 1 d D ! ! 2 E ð1Þ jrðtÞ -rð0Þ j lim D ¼ 2Nd t!¥ dt ! where rðtÞ is the position of a tagged species at time t and the angular brackets indicate an ensemble average. The parameter Nd is the dimensionality. Self-diffusivity of water, computed as an average over all the simulations conducted in this study was 4.05 ( 0.26 10-9 m2/s. The corresponding self-diffusivity of toluene was found to be 2.7 ( 0.3 10-9 m2/s. Although the diffusivity of water is in reasonable agreement with reported values in literature, the toluene diffusivity obtained in these binary film-like systems was
Figure 6. (a) Schematic depiction of the calculation of angle of inclination R between vector a (in the aromatic plane of asphaltene molecule) and vector n (normal to the toluene-water interface in xdirection). (b) Histogram of the relative angular distribution of four types of asphaltene molecules (number fraction, κi =Ni/Ntotal) represented as sin(R).
3.4. Orientation of Polyaromatic Rings with Respect to Toluene-Water Interfaces. The orientation of the asphaltene aromatic planes relative to the toluene-water interface is studied here. Simulations 1-4 all have identical number of asphaltenes, water, and toluene, which corresponds to the same point in the phase space (N, p, κi), where N is the number of molecules, p is the pressure, and κi is the composition (number fraction of species). Figure 6a shows the methodology of calculating the inclination angle R formed by the vector aB, lying on the aromatic plane of the asphaltene molecule and vector nB, which is normal to the toluene-water interface. Such a definition of the angle R allows the interpretation of the orientation of the polyaromatic regions of the asphaltene molecules with respect to the interface. The direction of the x-axis is normal to the interface. To obtain the angle of inclination, we calculated the average value of sin(R) and its standard deviation for each type of asphaltene molecule employing the simulation trajectories. The averaging was performed over the last nanosecond for each simulation. When sin(R) f 1, the aromatic planes of the asphaltenes are perpendicular to the toluene-water interface (or parallel to
(33) Zhang, L.; Greenfield, M. L. Energy Fuels 2007, 21, 1102–1111.
5033
Energy Fuels 2009, 23, 5027–5035
: DOI:10.1021/ef9004576
Kuznicki et al.
(circles) mean square displacements of the model continental asphaltenes in different binary solvent mixtures. As a reference, the mean square displacements of the C asphaltenes in pure toluene11 are shown in Figure 7a. Figure 7b and 7c show the mean square displacements of the C (simulation 1) and the AC (simulation 2) models, respectively, in the toluene film. One can represent the slopes of these transverse and lateral components of the mean square displacements versus time plots as transverse and lateral diffusivities, respectively. The representative straight lines (dashed red lines) are depicted in these figures to highlight the slopes of the lateral and transverse displacements. In case of asphaltenes in pure toluene (Figure 7a), the diffusion is relatively isotropic (slopes of lateral and transverse displacement plots against time are nearly identical). However, when the asphaltenes are trapped in the toluene film (Figure 7b and 7c), their long-term diffusion behavior appears to be dominated by lateral diffusion. In contrast, the transverse diffusion of the asphaltenes in the film is considerably suppressed. The difference between the transverse and lateral diffusivites is even more pronounced for the AC asphaltenes (Figure 7c). This type of anisotropic dynamics is generally observed in lipid bilayers where the lateral diffusion of the lipids along one monolayer (leaflet) is more dominant compared to transverse “flip-flop” diffusion.36,37 Furthermore, in simulations involving thin bilayers, the calculated diffusivities have been observed to be larger than experimental diffusivities, and such discrepancies could be attributed to several factors, including empirical force fields, limited box dimensions, and limited sampling.37 The anisotropy of the diffusion is strongly dependent on the asphaltene concentration in the toluene film. At lower concentrations of asphaltenes, the diffusivity becomes lower. This is evident by comparing the diffusivity of the VO-79 in simulations 4 and 5 (Table 1). In simulation 5, where the asphaltene concentration is five times lower than that of simulation 4, the diffusion coefficient is 3.8 ( 2.0 10-10 m2/s. 3.6. Discussion. One might find some similarities in the molecular architecture of VO-79 and some lipids (for instance, dipalmitoylphosphatidylcholine or DPPC).38,39 Both molecules have two long aliphatic tails suspended from a “head” region. For VO-79, the head region will be a polyaromatic ring, whereas for the DPPC lipids, these consist of the amine and phosphate groups. It is well-known that DPPC partitions at oil-water interfaces, with their ionic head regions embedded in the aqueous phase. With a similar molecular architecture, and the high polarity of the polyaromatic head of VO-79, one would perhaps expect these asphaltenes to exhibit similar properties as well. However, this is not the case in our simulations. The difference between VO-79 and lipids/surfactants is that the former is neutral and uncharged, whereas lipids have charged head groups. It seems that presence of charge (even in zwitterionic systems) is essential in creating a surfactant-like activity necessary for structuring of lipid head groups at aqueous interfaces. A similar behavior is evident in case of the model asphaltenes
Table 2. Fitted Parameters of a Normal and a Log-Normal Distribution Curve to the Simulation Based Distribution of sin(R)a normal distribution log-normal distribution simulation asphalt model mean 1 2 3 4
C AC TAC VO-79
0.81 0.71 0.70 0.80
std. error
mean
variance
0.15 0.14 0.20 0.08
0.813 0.712 0.707 0.804
0.027 0.032 0.058 0.006
a Where R is the angle between the aromatic plane of the asphaltenes and the unit surface normal to the toluene-water interface (see Figure 6).
about three times larger than the corresponding diffusion coefficients reported in literature.11,34,35 It is interesting to note that when we conducted simulations using pure toluene in our previous study,11 its diffusion coefficient was obtained as 0.89 10-9 m2/s, a value that is in excellent agreement with toluene diffusivity reported in literature. This indicates that thin films of toluene trapped within water phases exhibit a greater diffusion compared to the toluene molecules in an isotropic cubic simulation box. This may not exclusively be an artifact of the simulations since the molecular models, topologies, and other simulation parameters were unaltered in these studies. In our previous study, we calculated the diffusion coefficient of the model C asphaltenes (the same model as in Figure 1a) in pure solvents (water and toluene).11 The reported diffusivity of the C model asphaltenes in toluene was equal to 1.2 ( 0.8 10-10 m2/s. This is in good agreement with the reported diffusion coefficient of asphaltenes obtained by nuclear magnetic resonance (NMR).7 This agreement is encouraging since the concentration of the asphaltenes used in our earlier simulations was in the range of 0.34 to 1.2 g/L (above the CNAC of the asphaltenes), and the asphaltene MW was in the vicinity of 700 (same MW range used by Lisitza et al.7). These asphaltene diffusivities, although subject to considerable uncertainties owing to the small number of molecules (24) employed in the accumulation of the averages, provide some confidence in the ability of these MD simulations to provide reasonable predictions of the short-term dynamics of these model molecules. In the present study, however, the diffusion of the asphaltene molecules occurs within the toluene film. Therefore, the diffusion becomes considerably anisotropic in this thin film. It is noted that the diffusion coefficients calculated for these model asphaltenes (Table 1) are somewhat larger than the corresponding diffusivities in bulk (single phase) solvents.7,11 For instance, in simulation 4, the diffusivity of VO-79 is 4.7 ( 1.5 10-10 m2/s. The corresponding bulk diffusivity of an asphaltene that is approximately the same size as VO-79 was reported by Lisitza et al. as 2.85 10-10 m2/s.7 Clearly, as noted in case of toluene, we once again observe an increase in diffusivity when the molecules are trapped in a thin aromatic film between water layers. An interesting insight into the molecular dynamics of the asphaltenes in the toluene films is gained when one resolves the mean square displacement of these asphaltenes into two components, one normal to the toluene-water interface (the x-displacement), and the lateral displacement. Figure 7 depicts the total (solid line), lateral (triangles), and transverse
(36) Falck, E.; Rog, T.; Karttunen, M.; Vattulainen, I. J. Am. Chem. Soc. 2008, 130, 44–45. (37) Flenner, E.; Das, J.; Rheinstadter, M. C.; Kosztin, I. Phys. Rev. E 2009, 79, 011907. (38) Soliman, W.; Bhattacharjee, S.; Kaur, K. Biochim. Biophys. Acta-Proteins Proteomics 2007, 1774, 1002–1013. (39) Elmore, D. FEBS Lett. 2006, 580, 144–148.
(34) Bhide, S.; Berkowitz, M. J. Chem. Phys. 2005, 123, 224702. (35) Lee, S.; Lee, H.; Pak, H. Bull. Korean Chem. Soc. 1997, 18, 478– 484.
5034
Energy Fuels 2009, 23, 5027–5035
: DOI:10.1021/ef9004576
Kuznicki et al.
and align with respect to the toluene-water interface. In addition to π-stacking (responsible for face-to-face orientation of model asphaltenes within an aggregate), other types of interactions, such as the formation of hydrogen bonding at the toluene-water interface, takes place as a result of various (potentially charged) heteroatoms present in these asphaltene molecules. According to Xu et al.,40 interactions that explicitly involve localized polar functional groups (such as hydrogen bonding) are likely the most important when the binding of water to asphaltenes is considered. This formation of hydrogen bonds overpower π-stacking and results in location of AC and TAC asphaltenes at the interface whereas C and VO-79 do not exhibit any interfacial activity. 4. Concluding Remarks The MD simulations were conducted to elucidate the aggregation and partitioning of model asphaltenes in a binary solvent system, and relate these behaviors to their molecular structures. The study specifically aims to elucidate how the competition and synergy between asphaltene self-aggregation and interaction with toluene-water interfaces influence the aggregate morphology, position, and orientation relative to the toluene-water interfaces. Our simulations show that the model asphaltenes studied (C, VO-79, AC, and TAC), all tend to form aggregates above a critical concentration, which resemble the asphaltene nanoaggregates reported in several earlier studies. A key structural feature of such aggregates is the presence of stacked polyaromatic rings. The overall dynamics of the simulated systems involved a rapid phase separation of the solvents;toluene and water; into two film-like phases, followed by partitioning of the asphaltenes into the toluene phase. The long-term evolution of the aggregate nanostructures and their positions in the toluene films were different depending on the presence or absence of charge on these asphaltene models. Continental and VO-79 asphaltenes were preferentially positioned in the center of the toluene film (core-bound), showing no interfacial activity. In contrast, the AC and TAC molecules were interface-bound, where the charged terminal groups (COO-) of these molecules seemed to form hydrogen bonds with the water molecules at the toluene-water interface. It was observed that the polyaromatic ring regions of all asphaltene molecules prefer to position themselves perpendicular to the toluene-water interface. Finally, the simulations depict a highly anisotropic diffusion of the asphaltenes in the toluene films, particularly at high concentrations. The dominance of lateral diffusion of the asphaltenes seem to mimic lipid diffusion characteristics encountered in bilayers and membranes.
Figure 7. Variation of the mean square displacements of asphaltenes in different solvent systems. The straight dashed lines highlight the slopes of the long-term lateral (triangles) and transverse (circles) mean square displacements. When the two slopes are identical, the lateral and transverse diffusivities are similar. (a) Isotropic diffusion of C asphaltenes in a cubic periodic box of toluene. (b) Diffusion of Casphaltenes in toluene film (simulation 1). (c) Diffusion of AC asphaltene (simulation 2) in toluene film.
Acknowledgment. The authors gratefully acknowledge financial support from the Canada Research Chairs (CRC) program and Natural Sciences and Engineering Research Council of Canada (NSERC).
studied here, where absence of charge sites in the polyaromatic rings render these “head regions” surface inactive. The asphaltene-toluene and asphaltene-water pair correlation functions and orientations indicate that different chemistries and polarities in the four types of model asphaltenes studied can affect the ways these molecules aggregate
(40) Xu, Y.; Dabros, T.; Hamza, H.; Shefantook, W. Pet. Sci. Technol. 1999, 17, 1051–1070.
5035
