Understanding the Aggregation of Model Island and Archipelago Asphaltene Molecules near Kaolinite Surfaces using Molecular Dynamics

The solubility of asphaltenes in hydrocarbons changes with pressure, composition, and temperature, leading to precipitation and deposition, thereby causing one of the crucial problems that negatively affects oil production, transportation, and processing. Because, in some circumstances, it might be advantageous to promote asphaltene agglomeration into small colloidal particles, molecular dynamics simulations were conducted here to understand the impacts of a chemical additive inspired by cyclohexane on the mechanism of aggregation of model island and archipelago asphaltene molecules in toluene. We compared the results in the presence and absence of a kaolinite surface at 300 and 400 K. Cluster size analyses, radial distribution functions, angles between asphaltenes, radius of gyration, and entropic and energetic calculations were used to provide insights on the behavior of these systems. The results show that the hypothetical additive inspired by cyclohexane promoted the aggregation of both asphaltenes. Structural differences were observed among the aggregates obtained in our simulations. These differences are attributed to the number of aromatic cores and side chains on the asphaltene molecules as well as to that of heteroatoms. For the island structure, aggregation in the bulk phase was less pronounced than that in the proximity of the kaolinite surface, whereas the opposite was observed for the archipelago structure. In both cases, the additive promoted stacking of asphaltenes, yielding more compact aggregates. The results provided insights into the complex nature of asphaltene aggregation, although computational approaches that can access longer time and larger size scales should be chosen for quantifying emergent meso- and macroscale properties of systems containing asphaltenes in larger numbers than those that can currently be sampled via atomistic simulations.


INTRODUCTION
Asphaltenes, the heaviest and most polar fraction of crude oil, pose a flow assurance problem in the oil industry due to their tendency to form aggregates upon changes in temperature, pressure, and the physical-chemical environment.−5 Asphaltenes are generally defined as a solubility class depending on their solubility in aromatic solvents (such as toluene and xylene) and insolubility in aliphatic solvents (such as n-pentane and nheptane). 6,7As a consequence, their molecular structure is highly variable, and it strongly depends on the crude oil of interest.
In general, asphaltenes are represented either by the island structure, consisting of polycyclic aromatic hydrocarbon (PAH) rings attached to aliphatic chains, or by the archipelago structure, which consists of several aromatic cores loosely bonded by aliphatic chains.Chemical analysis has shown that PAH cores can contain heteroatoms such as nitrogen, oxygen, and sulfur, which often control aggregation, 8−11 and sometimes lead to hydrogen bonds. 9,12o interrogate the mechanism of asphaltene precipitation, molecular simulation studies have been attempted.Of these, several molecular dynamics (MD) studies showed that asphaltenes aggregate assuming face-to-face (parallel), Tshaped, or offset stacking structures.Interactions between PAH cores (π−π stacking) are frequently identified as the main driving force for aggregation.The size of the PAH cores controls the formation of π stacks and other aggregate types.−17 In addition to aggregation in the bulk, many practical applications are also affected by asphaltene adsorption on surfaces, a process likely controlled by both surface features and asphaltene characteristics.Adsorption of asphaltenes has been documented on mineral surfaces such as calcite, mica, and silica, 18−26 which affects both asphaltene aggregation and wettability of the mineral surfaces.Liu et al. 20,27 concluded from their experiments and MD simulations that polar silica surfaces are more attractive to asphaltene than nonpolar silica.They showed that the arrangement of heteroatoms on asphaltene contributes to the strength and mechanism of their adsorption.Xiong et al. 28 determined that the asphaltene model they used adsorbs on silica surface in heptane through hydrogen bonds with the hydroxyl groups on the surface, whereas in toluene, aggregates form and mostly stay in the bulk oil phase.Compared to silica, adsorption on muscovite is weaker and the molecules adsorb somewhat perpendicularly to the surface by interacting with the polar groups in the side chains or the edge of the PAHs. 29Mohammed and Gadikota 18,30 revealed that asphaltenes adsorb as small nanoaggregates on calcite and silica surfaces compared to the aggregates observed in bulk systems.The authors also demonstrated that CO 2 increases the aggregation of asphaltenes while displacing them from calcite due to the strong affinity between CO 2 and the surface.Fang et al. 31 observed, via MD simulations, that CO 2 causes asphaltenes to first form dimers and then aggregate and adsorb onto silica surfaces.Increasing pressure was found to weaken the interactions between asphaltenes and their adsorption on silica.Li et al. 32 experimentally observed CO 2 -induced asphaltene deposition on quartz.The aggregation and deposition were enhanced by increases in the temperature, resulting in a less water-wet surface.Li et al. 33 also reported that flooding with CO 2 (as opposed to flooding with methane and propane) leads to high asphaltene aggregation and found that the degree of aggregation increases with decreasing temperature.On quartz, the time taken for Violanthrone-79 (VO-79), a PAH frequently used as a model for asphaltene, 34 to adsorb on the surface increased with temperature, even though the number of stably adsorbed VO-79 molecules remained relatively constant. 35Moreover, the adsorption was more stable in heptane than in toluene, which allowed the PAH to desorb back into the bulk. 35Nassar et al. 24 found that asphaltene adsorbs more strongly on acidic alumina surfaces than on basic or neutral ones.These studies highlight the important contribution of the chemical composition of mineral surfaces in determining asphaltene precipitation.
To control and possibly prevent asphaltene precipitation, chemical inhibitors and dispersants are commonly used.Simulations have been utilized in identifying the mechanistic behavior of these additives.For example, employing MD simulations, Headen and Boek 36 showed that the presence of 50 wt % limonene reduces the aggregation of asphaltene in CO 2 .−40 According to Aminzadeh et al., 37 nonylphenols slightly reduce the size of asphaltene aggregates, with stronger effects resulting from nonylphenol added at the nucleation stage as opposed to being added earlier.Further increases in the concentration of nonylphenol can lead to their aggregation, a result ascribed to high dipole moments and the ability to form hydrogen bonds.Octylphenol also reduces the interaction energies and hydrogen bonds between asphaltenes.However, in terms of aggregation number, this inhibitor only delays asphaltene aggregation, rather than preventing it. 38The experiments reported by Lin et al. 40 showed that alkylphenols with longer alkyl tails increase the steric repulsion between aggregates, although the dispersants reduce the aggregate size leading to higher deposition rate.Strong repulsive forces between asphaltenes lead to the formation of "softer" aggregates that are readily eroded when shear flow is applied after deposition.Using MD, Jiang et al. 41 found that dodecylbenzenesulfonic acid (DBSA) reduces the rate and extent of asphaltene aggregation by reducing π−π stacking.When added after aggregation, DBSA weakens the aggregates by breaking the polar group bonds between the molecules but not the π−π stacking.Another study explored the effects of temperature on the performance of chemical inhibitors dodecyl benzenesulfonic acid (DBSA), benzoic acid (BA), and salicylic acid (SA) in the temperature range between 25 and 65 °C.DBSA was found to be the most effective inhibitor due to its acid−base interactions with asphaltenes.The increase in temperature accelerated asphaltene aggregation, which reduced the effectiveness of inhibitors.Salicylic acid caused the aggregates to increase at 65 °C as a result of strong hydrogen bonds. 42The results obtained for salicylic acid are comparable to those reported by Madhi et al., 43 where increasing the concentrations of cetyltrimethylammonium bromide (CTAB) and sodium dodecyl sulfate (SDS) led to self-assembly of the inhibitors, reducing their effectiveness.These studies, in general, show that the efficacy of a chemical inhibitor is determined by the strength of the interactions between asphaltenes and inhibitors relative to those between the two inhibitors.How this balance will be affected by the presence of a surface has not yet been fully quantified.
Despite this extensive body of work, the mechanism of action of asphaltene inhibitors is not well understood.For example, Kuang et al. 44 suggested that although conventional dispersants reduce the size of asphaltene aggregates, they do not necessarily alleviate asphaltene deposition.This leads to the possibility that enhancing the aggregation of asphaltenes into small structures could be a better strategy for preventing large plugs.Accordingly, the research presented here explores, on a molecular level, the use of a hydrocarbon chain to increase the aggregation of two asphaltene structures (island and archipelago), in combination with the proximity to a kaolinite surface.Kaolinite, a clay, is chosen because it is one of the common minerals found in reservoirs. 45,46Kaolinite is a layered aluminosilicate with the chemical formula Al 2 O 3 2SiO 2 • 2H 2 O.It has a 1:1 uncharged layered structure, consisting of alternating sheets of silica (SiO 4 ) tetrahedra and octahedral alumina oxyhydroxides, joined by apical oxygen atoms.−49 We employ MD simulations at the atomic resolution for this investigation.
The outline of this article is as follows: In the next section, we show the hydrocarbon models used and describe the force fields, initial configurations, and simulation algorithms.In section 3, we analyze the results by quantifying properties such as aggregate sizes, radial distribution functions, angles between the asphaltene molecules, and entropy analysis.Finally, we summarize in the conclusions our main findings and suggest improvements to our method.
Energy & Fuels 2. SIMULATION METHODOLOGY Molecular Models.In Figure 1, we present the asphaltene models used in this study.Violanthrone-79 (712.9 g/mol), 29,50 which consists of a large central poly aromatic core connected to two aliphatic side chains, represents the island structure (indicated as ASPH in what follows, shown in Figure 1a).This compound has been used experimentally to conduct controlled tests; it shows a tendency to aggregate in heptane to a larger extent than in toluene, which supports the choice of Violanthrone-79 as a model for asphaltene. 51As for the representative archipelago (ARCH) structure, a molecule with two aromatic cores, with molecular weight 792.2 g/mol, was chosen (Figure 1b). 52Although there have been several models proposed for asphaltenes, we chose to focus on one island model and one archipelago model to test the effects of the cyclohexane oligomer as a proof of principle for the potential of a putative chemical additive to increase the aggregation of asphaltenes with different structures.Of note, Violanthrone-79 is a real polycyclic aromatic compound, allowing our results to potentially be compared to experiments.Because previous research showed that asphaltenes are almost insoluble in cyclohexane, 53 we constructed a cyclohexane oligomer (Figure 1c) to investigate its effects on the aggregation of the asphaltene models.
Although crude oil is a mixture of several hydrocarbon components, for the purpose of this study, the oil model was represented by only toluene.By definition, toluene is a good solvent for asphaltenes.It is chosen here because the aim is to test whether additives can promote asphaltene aggregates formation.In other words, we assume this system represents a starting point to quantify the increased aggregation of asphaltenes due to chemical additives.
Simulation Systems.We constructed three simulation systems.Two of these (Systems 1 and 4) contain only toluene and asphaltene to benchmark the other simulations.The second group (Systems 2 and 5) and third group (Systems 3 and 6) of systems, as shown in Table 1, contain varying amounts of toluene and cyclohexane oligomer.In all systems, we used 12 asphaltene molecules to maintain the percentage by mass of asphaltene at ∼8 wt %, 13,16,54,55 which is representative of asphaltene fraction in crude oil. 56,57Because of the difference in the molecular masses of the island and archipelago asphaltenes, different numbers of molecules of toluene and cyclohexane chain were inserted in our simulated systems.For example, since the archipelago model used has a slightly higher molecular mass, it requires more toluene and cyclohexane chains to maintain the three systems with comparable weight fractions.
For all systems shown in Table 1, we conducted simulations in the bulk, and near a kaolinite surface.In the latter case, we constructed the simulation box by placing a thick slab of asphaltene, toluene, and cyclohexane oligomers on the kaolinite surface with dimensions 5.198 × 8.982 × 1.410 nm 3 parallel to the X−Y plane, resulting in a simulation box of size 5.198 × 8.982 × 9.500 nm 3 .The asphaltene molecules were inserted randomly into the box, ensuring there is little aggregation at the start of the simulation.Figure 2 shows the snapshot of the initial configuration of the system.The hydrocarbon phase is sandwiched between two parallel solid surfaces, at a distance ranging between 3.7 and 4.5 nm.For simple fluids, confinement effects are pronounced within one or two molecular layers from the solid−liquid interface.The size of the system is chosen to be larger than that, but we cannot exclude a priori that confinement impacts the results.Unfortunately, at the resolution considered for this study, larger distances between the kaolinite surfaces quickly become prohibitive because of increased computational costs.The   geometry chosen prevents the formation of a liquid−gas interface, which could compete for asphaltene adsorption compared to the solid−liquid interface.
The bulk simulation boxes were constructed in a similar manner, except without the kaolinite surface.The size of the simulation box used for bulk simulations was 5.198 × 8.982 × 6.500 nm 3 .The resultant density of the hydrocarbon was approximately 3.39 molecules/nm 3 .In all systems, periodic boundary conditions were applied in all directions in the simulation box.
Force Fields.The Optimized Potentials for Liquid Simulations All Atom (OPLS-AA) force field 58 was used to describe all hydrocarbons considered, following numerous related works in the literature. 13,14,56,59,60OPLS-AA has been shown to reproduce reasonably well experimental data for aromatic liquids. 61,62The hydroxylated kaolinite surfaces were modeled using the CLAYFF force field, 63 following prior research on clay minerals from our group. 64,65−68 Prior studies confirmed that by combining CLAYFF and OPLS-AA force fields it is possible to achieve simulation results consistent with experimental observations.For example, Schampera et al. 69 studied the interactions of organic cations on a montmorillonite clay surface, and their results were consistent with experimental XPS data.As another example, Ren and Liu 70 and Phan et al. 71 studied the structure of ethanol on alumina surfaces and found that the combination of these force fields compares well with experimental data.
Dispersive forces were modeled by the 12−6 Lennard-Jones (LJ) potential.Electrostatic forces were described by the Coulomb potential.The LJ parameters for unlike atomic interactions were obtained by using the geometric mixing rules, commonly used with the OPLS family of force fields.In all nonbonded interactions, the cut off radius for short-range interactions was set to 14 Å.Corrections for long-range electrostatic interactions were obtained by the particle mesh Ewald (PME) algorithm. 72,73lgorithms.Atomistic MD simulations were conducted using the GROMACS package, version 5.1.4. 74,75After energy minimization, we carried out NVT canonical simulations for 1 ns to relax the initial configuration with the positions of the kaolinite surface kept fixed.The simulations were then conducted in an isobaric−isothermal (NPT) ensemble for 200 ns to ensure the densities and energies reached a stable state and enough for aggregation to be observed within a reasonable time scale.A time step of 1 fs was used for all simulations.During the NPT simulations, the kaolinite surface was restrained with a force constant of 1000 kJ/mol nm 2 while asphaltenes, toluene, and the cyclohexane oligomer were allowed to move.
Subsequently, the time-averaged properties of the system (such as densities and radial distribution functions) were analyzed by conducting two additional 4 ns simulations for each system.The coordinates of the simulated systems were extracted every 0.5 ps to construct the simulation trajectories.The results were then quantified using the techniques described in the Results section.We set the temperature to either 300 or 400 K, with a relaxation time of 100 fs.The temperature was controlled by a Berendsen thermostat.−78 To probe the possible mechanisms responsible for agglomeration, only Systems 1 and 3 (Table 1) were simulated at 400 K.The pressure was maintained at 15 MPa by the Berendsen barostat.This pressure represents that found in reservoirs. 36,79In the systems with a kaolinite surface, pressure coupling was applied only along the Z direction of the simulation box, perpendicular to the surface.This ensures that the surface area (X and Y dimensions) is constant for all systems.

RESULTS
Density Profiles.The density profiles of the center of mass of toluene, island asphaltene (ASPH), and the cyclohexane oligomer (HEX), computed along the perpendicular (Z) direction at 300 and 400 K, are presented in Figure 3.The top plane of the hydroxyl ions represents the reference position for the vertical distance Z. From the results, we observe that at both temperatures, toluene adsorbs closest to the kaolinite surface, forming three layers with peaks at distances corresponding to 2.1, 3.7, and 7.3 Å from the surface.At 300 K, the first toluene layer achieves almost twice the density of the second and third layers.The island asphaltene adsorbs further from the surface from ∼3.1 Å up to 12 Å, with a pronounced peak at 3.9 Å, suggesting that toluene forms a "protective" layer on kaolinite.We can infer that toluene has stronger interactions with the kaolinite surface than the island asphaltene does.In contrast, the cyclohexane oligomer (HEX) yields no obvious layering on kaolinite.Our simulations show that adding HEX has no distinct effect on the structure of    Energy & Fuels adsorption of the poly aromatic hydrocarbon on the mineral surface. 35lthough the center of mass of the asphaltenes seems further from the kaolinite surface than toluene, we expect that the presence of oxygen atoms could lead to bonds with the hydroxyl groups on the surface.To confirm this, we calculated the atomic density distribution of 2 oxygen atoms present in the asphaltene structure.The results are highlighted in Figure 4. We observe a sharp peak at ∼1.5−1.7 Å, which is closer to the surface than the density peaks computed for the centers of mass of both asphaltene and toluene.The height of the peak increases with temperature, in agreement with the observations from the center of mass density profiles of asphaltene.From the positions of these oxygen atoms, we infer that the main interactions of asphaltene with the surface are due to the presence of heteroatoms, which could potentially form hydrogen bonds with the hydroxyl ions on the mineral surface.
In Figure 5, we provide the center of mass density profiles obtained for the systems containing the archipelago-structured asphaltene (ARCH).The positions of the toluene and HEX peaks correspond approximately to the positions shown in Figure 3.However, the archipelago-structured asphaltene adsorbs somewhat closer to kaolinite, with a distinct peak at 2.5 Å compared to the island asphaltene, whose first peak appears at 3.9 Å.A possible explanation is that the presence of sulfur, oxygen, and nitrogen atoms in the archipelago structure increases electrostatic and van der Waals interactions with the kaolinite surface (note in Figure 1 that the island model used here does not contain sulfur nor nitrogen atoms).At 300 K, a second asphaltene peak can be observed at 12.9 Å (System 4) and 12.1 Å (System 5); this peak disappears upon adding HEX or upon increasing the temperature to 400 K.This is due to the asphaltene packing at a closer distance to the surface, as observed from the noticeable increase in the height of the first asphaltene peak at 2.5 Å.It is worth noting that at 400 K, HEX does not significantly affect the density distribution of asphaltene near the surface.
To probe the mechanisms responsible for the preferential adsorption, we computed the density profiles of the nitrogen and oxygen atoms found in the archipelago asphaltene structure.The results are presented in Figure 6, where it is shown that the relevant density peaks are found between 1.5 and 1.7 Å from kaolinite.These results highlight the interactions between the heteroatoms and clay mineral surfaces.It is worth noting that we did not observe a strong influence of HEX on the positions or intensity of the peaks, which could be due to the complex nature of the aggregation of archipelago asphaltenes.
Cluster Size Analysis.To quantify the aggregation of asphaltene molecules, we estimated the average cluster size by calculating the number of asphaltene molecules in each aggregate and taking the average of the results obtained over the simulated trajectories.We also quantified the maximum number of asphaltenes in an aggregate.Following several studies in the literature, 13,18,30,38 asphaltene molecules are considered clustered if the distance between them (distance of closest approach) is less than 3.5 Å.As shown in Figure 7a, at 300 K, the average cluster size of the island asphaltene (ASPH) increases as the number of HEX oligomers in the system increases.In the last 50 ns of the simulation, the average aggregate size in the system with only toluene (System 1) is ∼6−7.This increased to ∼7−8 and ∼8−9 in Systems 2 and 3, respectively.In System 3, asphaltenes yield one unstable aggregate, which breaks off repeatedly.The results in Figure S1, panel a, in the Supporting Information show that adding   1.
HEX leads to an increase in the maximum number of asphaltenes seen in any aggregate.
When the temperature is increased to 400 K, Figure 7b shows that there are more frequent fluctuations in the aggregate size, which is likely a consequence of the increase in the kinetic energy within the system.It is also interesting to note that the average size of the aggregates in systems 1 and 3 are ∼6.3 and ∼7.1, respectively.These results suggest that increasing the temperature increases the number of HEX required to cause a marked increase in the average cluster size, a result which is likely due to an increase in disorder of the system.
We observe trends similar to those just discussed in the simulated bulk systems.In Figure 8, we present the cluster size of ASPH asphaltenes in the bulk.In System 1 (only toluene), the average size of the aggregates is ∼3−4, which increases to 4−5 with the addition of 22 HEX oligomers (System 2) and 5.8 with the presence of 87 HEX oligomers (System 3).In System 1, the maximum cluster size stabilizes at 6 asphaltenes, while this value increases to 7 and 8 in Systems 2 and 3, respectively.These results suggest that asphaltene aggregates in bulk systems are more stable and, on average, smaller than aggregates near kaolinite.
In Figure 9, we provide cluster size analyses for archipelago asphaltene (ARCH) on kaolinite at 300 K.The effects due to the cyclohexane oligomer are similar to those discussed above for island asphaltene.In the last 50 ns, the average aggregate size increases from 6 to 10 when 100 HEX oligomers (System 6) replace some toluene in the system.Correspondingly, the maximum aggregate size increases from ∼7.9 to ∼10.9.The aggregates observed here are slightly larger than the aggregates formed by the island asphaltene.
In the bulk, the clustering results show that all the asphaltenes aggregate (Figure S2), even without HEX.The same is observed when the temperature is increased from 300 to 400 K in proximity of kaolinite.Snapshots of the fully aggregated asphaltenes are reported in the Supporting Information (Figure S3).We attribute this result to the difference in the size of the poly aromatic cores, length of the side aliphatic chains, and the presence of heteroatoms, compared to the relevant features of the island asphaltene structure.In the island asphaltene, the aromatic core contains 8 aromatic rings, whereas the archipelago asphaltene has 2 poly aromatic cores each of 5 and 6 aromatic rings, which could increase the π−π interactions.
Radial Distribution Functions (RDF).In Figure 10, we report the RDF between the center of mass of the aromatic cores of the island asphaltene (ASPH) at 300 and 400 K for the systems in proximity of the kaolinite surface.Because only layers occupied by asphaltenes (up to ∼1.5 nm), shown by the atomic density distributions, were considered for the RDF analyses, the systems are nearly planar.Because interfacial systems are not isotropic, the RDFs were computed in 2 dimensions.In system 1, panel a, the first most prominent peak is found at a radial distance of ∼5.3 Å.This suggests offset stacking of the aromatic cores or a slant orientation between asphaltene pairs.Carauta et al. 80 reported from experiments that asphaltene dimers bind face−face at a distance of ∼3.6 Å in heptane and ∼5 Å in toluene.The peak distance is slightly greater than the π−π stacking distance of 5 Å.In the presence of HEX, the height of the first peak increases and is observed at a slightly closer radial distance of 4.8 Å, indicating stronger π−π interactions and a greater tendency to form aggregates.For all the systems, there are additional peaks observed at 10.6  The RDF at 400 K (panel b) shows a similar trend due to the HEX increasing the interactions between asphaltenes.In System 1, the most pronounced peak is found at a radial distance of ∼6.3 Å, which agrees with the possibility that at a higher temperature, the island asphaltene forms loose aggregates.However, in the presence of 87 HEX, the peak shifts to a smaller radial distance of ∼4.8 Å even at the higher temperature simulated.If we compare the RDF of the same system at two temperatures, we find that increasing the temperature reduces the probability of observing asphaltene cores close to each other.
The toluene-ASPH RDF, presented in Figure 11, reveals the effects of HEX on the interactions between asphaltene and toluene.Interestingly, the results show that the presence of the HEX in the layers occupied by asphaltene has no strong effects on the structure of toluene around asphaltene at the two temperatures considered, even though in some systems, some HEX are present in the same layer as asphaltene and toluene.At both temperatures (panels a and b), the first and second peaks are located at ∼6.5 and 11.3 Å, respectively, although the height of the peak is slightly higher at the lower temperature (panel a).
For bulk systems containing the island asphaltene, Figure S4 of the Supporting Information demonstrates that the threedimensional RDF for asphaltene pairs is similar to the results obtained for the systems simulated in proximity of kaolinite; HEX is also found to increase the height of the peak located at 4.8 Å.Moreover, the toluene-asphaltene three-dimensional RDF in the bulk is unaffected by the presence of HEX.
For the archipelago asphaltene (ARCH), we computed RDFs between the center of mass of the asphaltene pairs.Because the multiple aromatic cores and the heteroatoms present in the archipelago-like asphaltenes lead to more complex stacking and aggregation, we did not use the RDF between the aromatic cores to analyze stacking between asphaltene molecules.The results, in Figure 12, show that the presence of the cyclohexane chain leads to closer interactions between the archipelago asphaltene.In Systems 5 and 6 at 300 K (panel a), we observe multiple peaks, with the first peak found at ∼1.9 Å, whereas in system 4 with only toluene, there are no sharp peaks.This suggests that compact aggregates form in the presence of the cyclohexane chain.A similar trend is found at 400 K (panel b), where the first broad peak of system 4 is found at 5.1 Å, and upon the replacement with the cyclohexane oligomer, the peak shifts to a smaller distance of 2.8 Å.The asphaltene−toluene RDF, presented in Figure S5 in the Supporting Information, shows no visible differences between the system with only toluene and the systems with toluene and cyclohexane oligomer.In the RDFs at the two temperatures, there are no significant spikes, indicating that no well-defined aggregate forms between toluene and asphaltene molecules within the systems considered here.

Orientation and Number of π−π Contacts of Island Asphaltene Molecules.
To gain insights on the structure of the asphaltene aggregates, we calculated the angle formed between the poly aromatic planes of two asphaltenes in an aggregate on the kaolinite surface (illustrated in Figure 13).The angle ranges from 0 to 90°and could be face−face or Tshaped.In the face−face stacking, the center of aromatic cores could be aligned (parallel) or slightly displaced (offset).The orientation gives more information about the interactions between the aromatic cores of the asphaltene, with the parallel orientation being the strongest.
In Figure 14, we report the angles between island asphaltene pairs at close contact (as identified by the first peaks in the RDFs of Figure 10) at 300 (panel a) and 400 K (panel b) for the kaolinite systems.As the number of HEX oligomers increases, the asphaltene pairs adopt a parallel orientation, suggesting a more compact aggregate structure established through π−π interactions.In System 1 at 300 K (panel a),   1.
there are two peaks at 5.7°and 23.7°.As the number of HEX increases, the probability of the angles between asphaltene pairs being less than 20°increases, with a pronounced peak at 5.7°in System 3.These results show that HEX oligomers shift the stacking between asphaltene from offset toward either parallel or offset structures.In the first RDF sphere of System 1 at 400 K (panel b), there is the highest probability of finding two asphaltene molecules oriented in a way to form angles between 40−50°, corresponding to a more T-shaped orientation.However, adding HEX shifts the peak to a smaller offset angle at 9.9°.The difference in the angles observed corroborates the findings that increasing temperature yields less compact aggregates.We observed similar trends in bulk systems (Figure S6).
We also quantified the relationship between the distance between asphaltene pairs and the type of stacking (Figure 15).We observed that as the distance between the asphaltenes increases, the orientation tends toward offset or T-shaped.In System 1, asphaltene molecules found at distances less than 7.1 Å are mostly offset, whereas those found at larger distances show a preference for T-shaped orientations.
To complement this analysis, we quantified the number of π−π contacts between aromatic cores of the island asphaltenes.A contact was considered established when the distance between the center of mass of the aromatic cores is less than ∼4.8 Å (the radial distance at which peaks are observed on the RDF).The results in Figure S7 in the Supporting Information show that the HEX oligomers increase the number of compact π−π contacts in the aggregated asphaltenes.In the presence of the kaolinite surface, at 300 K (panel a), the number of contacts in the absence of HEX fluctuates between 0 and 1; this value increases to 2−4 in System 3. At 400 K (panel b), the number of contacts decreases, in agreement with our previous observations that increasing the temperature leads to less compact aggregates.In panel c, we find that the number of π−π contacts is higher in the bulk phase than in the presence of the kaolinite slab.We see in panel c that the number of contacts in System 1 is between 0 and 2, while in System 3, the value increases to 6, indicating more rigid aggregates in bulk systems compared to the systems near kaolinite.
Compactness of Archipelago Asphaltenes.The asphaltenes of the archipelago type show more complex aggregation because of their extended molecular structure.To evaluate their structure, we quantified the distance between nitrogen (N) and sulfur (S) atoms in individual asphaltenes (see Figure 16).The results reported in Figure 17 (panel a) show that in the proximity of kaolinite at 300 K, although the asphaltene seems to completely unfold, the distance between the N and S atoms reduces slightly from ∼51 Å in System 4 to ∼43 Å in System 6, yielding slightly more compact molecules.Increasing the temperature to 400 K (Figure S8) reduces these values to ∼34 Å, which is consistent with high degree of aggregation, for both systems considered.In contrast to the effects of temperature on the island asphaltenes, increasing temperature makes the archipelago asphaltenes more compact.The corresponding results for bulk systems are presented in Figure 17 (panel b).HEX leads to a reduction in the distance, but it is worth noting that the effect is not very pronounced near kaolinite.
We also computed the radius of gyration (R g ) to support the analysis above.In the bulk phase (Figure S9), R g fluctuates within a wider range than in proximity of kaolinite (Figure S10).On kaolinite at 300 K (Figure S10a), R g reduces from ∼38 Å in System 4 to ∼31 Å in System 6 upon adding HEX oligomers.At 400 K (Figure S10b), R g remains at ∼35 Å in both systems.
Interaction Energies and Entropy.In Table 2, we report van der Waals (ΔE VDW ) and electrostatic interaction potentials (ΔE ELEC ) for ASPH molecules in the systems studied.The interactions between asphaltenes are dominated slightly (∼55%) by electrostatic interactions, due to the presence of heteroatoms, which increase molecular polarity.The presence of HEX oligomers has greater effect on the van der Waal interactions, which complements the observations of increase in π−π interactions.The interactions of the archipelagostructured asphaltenes vary in a similar manner and are shown in Table S1 in the SI.
The configurational entropy of asphaltene in toluene in the presence of HEX oligomers was estimated by applying the   1. method described by Schlitter. 81We determined the entropy of 3 asphaltenes with similar mass percentage to those used in our systems.Even though the uncertainty is large as a consequence of the small number of molecules considered for this calculation, the results in Table 3 suggest that the presence of HEX oligomers reduces the entropy of the asphaltenes, indicative of a more ordered structure.As expected, increasing the temperature leads to more disordered systems (higher entropy), which is consistent with results discussed above.

SUMMARY
We performed molecular dynamics simulations to study the aggregation of asphaltene, a major component of crude oils, in the proximity of a kaolinite surface and in the bulk.As a proof of concept, to determine the effects of a cyclohexane oligomer on model island and archipelago asphaltenes, we analyzed the aggregation of asphaltene in toluene, in the presence or absence of the oligomer, at 300 and 400 K temperatures, at the pressure of 15 MPa.The results suggest that the cyclohexane oligomer slightly increases the level of aggregation of asphaltene.In the island asphaltene, this was evident by the increase in the number of parallel π−π orientations of the aggregates, whereas in the archipelago structure, we observed more compact asphaltenes.
We found that the aggregation of the island model was more pronounced near kaolinite compared to the aggregation in the bulk phase, while the aggregation of the archipelago model was less pronounced.It is possible that the difference in the interactions between asphaltene molecules and the kaolinite surface, mediated with the formation of a dense toluene layer near the solid surface, contributes to this result.Our energetic and entropic analyses align with our observations.It should, however, be noted that confinement effects cannot be ruled out, given the relatively small distance between the two kaolinite surfaces confining the system investigated.
Overall, our study contributes to our understanding of the aggregation of asphaltenes near the kaolinite clay surface and the influence of the structure on aggregation behavior.For future work, larger systems with more asphaltene models should be investigated to monitor in greater detail the stages of aggregation that we could not implement via MD simulations due to high computational costs.A suitable alternative could be large-scale, coarse-grained molecular simulations.Experiments such as dynamic light scattering could also be conducted    1.
to directly investigate the aggregation of asphaltenes and compare the results to simulation results. 82ASSOCIATED CONTENT * sı Supporting Information The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.energyfuels.3c00504.
Maximum cluster size of island asphaltenes as a function of time; average cluster size as a function of time for the archipelago asphaltene on kaolinite at 400 K and in bulk systems at 300 K; RDFs between center mass of island asphaltene and center of mass of asphaltene and toluene in bulk systems at 300 K; RDFs between the center mass of archipelago asphaltene and the center of mass of toluene on kaolinite surface; angle between the poly aromatic planes of the island asphaltene pairs; number of π−π contacts between the aromatic cores of island asphaltene pairs on kaolinite; distance between sulfur and nitrogen atoms for the archipelago asphaltene on kaolinite at 400 K; radius of gyration of individual archipelago asphaltene in kaolinite and bulk systems; asphaltene-asphaltene interaction energies between archipelago asphaltenes; representative simulation snapshots (PDF) ■ System compositions are shown in Table S2.
Energy & Fuels

Figure 2 .
Figure 2. Representative snapshots of the initial configuration of the systems with hydrocarbons on the kaolinite surface.Cyclohexane chain = purple.Toluene = gray.

Figure 3 .
Figure3.Centre of mass profiles along the Z direction, vertical from the kaolinite surface, for island asphaltene (ASPH), toluene (Tol), and cyclohexane oligomer (HEX) at (a) 300 and (b) 400 K.For system composition, we refer to Table1.

Figure 4 .
Figure 4. Atomic density profiles of two oxygen atoms (circled) located on island asphaltene (ASPH) along the Z direction, vertical from the kaolinite surface.Results are for simulations conducted at 300 and 400 K.

Figure 5 .
Figure 5. Centre of mass density profiles along the Z direction, vertical from the kaolinite surface, for the archipelago asphaltene (ARCH), toluene, and cyclohexane chain (HEX) at (a) 300 and (b) 400 K. System compositions are shown in Table1.

Figure 6 .
Figure 6.Atomic density profiles of oxygen atoms (panel b) and nitrogen atoms (panel c) found on ARCH along the Z direction, perpendicular to the kaolinite surface.Simulations results are obtained at 300 and 400 K.In panel a, the oxygen atom is highlighted by the red circled sphere, while nitrogen is highlighted by the blue circled sphere.

Figure 7 .
Figure 7. Average cluster (aggregate) as a function of time for the ASPH on the kaolinite surface at (a) 300 and (b) 400 K. System compositions are shown in Table1.

Figure 8 .
Figure 8.(a) Average cluster (aggregate) size and (b) maximum cluster size as a function of time for ASPH in the bulk systems at 300 K. System compositions are shown in Table1.

Figure 9 .
Figure 9. Average cluster (aggregate) size and (b) maximum cluster size as a function of time for ARCH on kaolinite at 300 K. System compositions are shown in Table1.

Figure 10 .
Figure 10.RDFs between the center of the aromatic cores of the island asphaltene (ASPH) pairs found in proximity of the kaolinite surface at (a) 300 and (b) 400 K. System compositions are shown in Table1.

Figure 11 .
Figure 11.RDFs between the center mass of island asphaltene (ASPH) and the center of mass of toluene in proximity of kaolinite at (a) 300 and (b) 400 K. System compositions are shown in Table1.

Figure 12 .
Figure 12.RDFs between the center of the archipelago asphaltene (ARCH) pairs at (a) and (b) 300 K and (c) 400 K. System compositions are shown in Table1.

Figure 14 .
Figure 14.Angle between the poly aromatic planes of the island asphaltene (ASPH) pairs at (a) 300 and (b) 400 K on the kaolinite surface.System compositions are shown in Table1.

Figure 15 .
Figure 15.Angle between the poly aromatic planes of the ASPH pairs versus distance between them in (a) System 1 and (b) System 3 on kaolinite surface.System compositions are shown in Table1.

Figure 16 .
Figure 16.Illustration of the distance between sulfur and nitrogen atoms calculated for ARCH molecules.

Figure 17 .
Figure 17.Distance between sulfur and nitrogen atoms calculated for ARCH (a) on the kaolinite surface and (b) in the bulk at 300 K. System compositions are shown in Table1.

Table 1 .
Composition of the Bulk and Kaolinite Systems Used in Simulations Conducted at 300 and 400 K a a ASPH represents the island model asphaltene, while ARCH is the archipelago model.Only systems 1, 3, 4, and 6 were simulated at 400 K.

Table 2 .
Asphaltene−Asphaltene Interactions between Island Asphaltenes (ASPH) on the Kaolinite Surface and in Bulk Systems aOnly short-range contributions to potential energies are considered.System compositions are shown in Table1. a

Table 3 .
Configurational Entropy of 3 Asphaltenes in Bulk Systems at 300 and 400 K a