# Aggregation of Asphaltene Subfractions A1 and A2 in Different Solvents from the Perspective of Molecular Dynamics Simulations

Orlando Villegas, Germain Salvato Vallverdu, Brice Bouyssiere, Sócrates

Acevedo, Jimmy Castillo, Isabelle Baraille

► To cite this version:

Orlando Villegas, Germain Salvato Vallverdu, Brice Bouyssiere, Sócrates Acevedo, Jimmy Castillo, et al.

HAL Id: hal-03996687

https://univ-pau.hal.science/hal-03996687

Submitted on 20 Feb 2023

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Aggregation of asphaltene subfractions A1 and A2 in different solvents from the perspective of molecular dynamics simulations

Orlando Villegas,†,‡ Germain Salvato Vallverdu,*,‡,¶ Brice Bouyssiere,‡,¶ Sócrates Acevedo,† Jimmy Castillo,† and Isabelle Baraille‡,¶

†Facultad de Ciencias, Escueal de Química, Universidad Central de Venezuela, Caracas, 1041, Venezuela

‡Universite de Pau et des Pays de l'Adour, E2S UPPA, CNRS, IPREM, UMR 5254, Pau France

¶International Joint Laboratory C2MC: Complex Matrices Molecular Characterization, Total Research & Technology, Gonfreville, BP 27, F-76700 Harfleur, France

E-mail: germain.vallverdu@univ-pau.fr

## Abstract

The aggregation behavior of model molecules of asphaltene subfractions A1 and A2 dissolved in heptane, toluene and THF were studied using molecular dynamics simulations. The proposed asphaltene molecular models are based on previously studied structures with two new models, including a highly aromatic model, with a prominent island-type molecule and another molecule with prominent archipelago-type architecture. The aggregation mechanisms in toluene, THF and heptane solvents were studied. The results in heptane and toluene were consistent with the solubility of asphaltenes and their subfractions in these solvents. The size of the aggregates is well correlated with aromaticity. When considering THF, large aggregates are broken down into smaller aggregates. This could lead to the mixture of high, medium, and low molecular weight distribution bands usually observed when GPC analyses are conducted on asphaltene samples in THF. The size distribution extracted from the simulations shows a bimodal distribution with profiles similar to the size distribution profiles usually found in GPC analysis for asphaltene samples in THF. The

distributions of dipole moments of the aggregates against the number of molecules in aggregates were constructed both in THF and toluene and reveal that the dipole moment of the aggregates vanishes when the number of molecules increases due to a random structure of the aggregates. The contributions of different molecular interactions to the aggregation mechanism, such as π-stacking, van der Waals and hydrogen bonds, are described.
# 1 Introduction

Asphaltenes are a family of compounds defined experimentally by solubility, according to the SARA classification as the fraction insoluble in n-alkanes (usually n-heptane) and soluble in aromatic solvents (such as toluene).¹ Asphaltenes are the heaviest and most polar components of crude oil. They can precipitate at any point in production or transport as a result of changes in temperature, pressure, or composition.² This simple method of defining classes of compounds based on their solubility led to complications in determining their molecular weight, molecular structure, and chemical properties: consequently, practical recommendations to ease production and reduce undesirable behavior are difficult to define. In this way, several studies have been carried out to evaluate different models and suggest theories to determine the behavior of the asphaltene fraction in solution. In the family of asphaltene molecules, we usually distinguish molecules with one (continental or island type) or several (archipelago type) aromatic carbon nuclei, including heteroatoms (mainly N, S and O) with alkyl side chains of various lengths grafted onto these nuclei.³,⁴

The aggregation mechanism of asphaltenes has been a matter of much discussion over the years. The molecular characterization of asphaltenes in the condensed phase is the point of interest in most of the research published in this field. Investigation of asphaltenes in crude oil is complex in particular because this is a multiphasic system. Asphaltenes exhibit properties derived from the theory of colloidal solutions or colloidal models (or suspension models), and phase behavior originating from multicomponent mixtures or the solubility model. This shows the dual nature of asphaltenes: partly dissolved and partly colloidal.⁵ In recent years, attention has been focused on the colloidal model.⁶ It has been shown that asphaltenes form nanocolloids in crude oil which are very similar to the nanoaggregates they form in toluene.⁶ The molecular characterization of asphaltenes requires a detailed understanding of the colloidal hierarchical structure. The Yen-

Mullins model predicts that the asphaltenes form small and dense nanoaggregates of approximately 6 molecules, which are assembled to form clusters composed of approximately 8 nanoaggregates.⁷ The leading forces of this aggregation mechanism are π-stacking interactions between aromatic cores of asphaltenes associated with van der Waals and dipole interactions. The size of the nanoaggregates is controlled by steric repulsion due to the presence of the alkyl substituents (side chains) connected to the molecular core.⁸,⁹ Gray et al.¹⁰ proposed a supramolecular assembly model, which includes more diversity in the considered molecular interaction and provides a more complete view of the physical-chemical behavior of the aggregates and their complexity. In particular, it includes the porosity of the aggregates and their ability to trap compounds. The self-association of asphaltene molecules in solution is one of their canonical properties.¹¹ Understanding how the molecules interact with each other to form aggregates, and then how the aggregates arrange themselves into flocs or clusters, and ultimately precipitate, is very important in determining the stability of the asphaltene fraction in heavy oil processing.¹²

The study of asphaltene subfractions helps us to better understand the composition of this family of compounds. It has already been observed that the use of p-nitrophenol makes it possible to split the asphaltene fraction into two subfractions called A1 and A2 which exhibit differences in solubility.¹³,¹⁴ A2 shows the same solubility as the asphaltenes in toluene, while A1 is very poorly soluble. A recent hypothesis proposes that the mechanism for the aggregation of asphaltenes in toluene solutions proceeds primarily by the formation of small and dense aggregates of A1 forming nuclei that promote phase separation in the medium, while the intercalation of A2 molecules and their arrangement at the periphery of these nuclei help to stop the phase separation, thus forming a lipophilic colloid.¹³⁻¹⁵ This hypothesis is consistent both with Gray's supramolecular and with the Yen-Mullins models. Recently our group showed that asphaltenes form aggregates which are stable in various solvents such as toluene
and THF. Changing the solvent will impact the width of the hydrodynamics volumes distribution.¹⁶

Molecular dynamics simulations provide an essential tool to increase our knowledge of the behavior of the asphaltene fraction by considering, at an atomic scale, the intermolecular interactions in relation to the molecular structure. This is useful to understand the behavior of asphaltenes in solution at a nanometer scale and to link the results with experimental results on macroscopic samples. The main drawback of this approach is the requirement for the molecular structure of the system to be considered. This becomes very difficult for asphaltenes because they are constituted by a family of compounds; however the general approach consists of considering a subset of molecular models that are representative of their behavior in solution.¹⁷⁻²¹ Nevertheless, this approach allows the design of molecular models with known chemical characteristics anduncouples the key factors one seeks to investigate. We have already worked with molecular models that provide us with the possibility to explain experimental behaviors associated with the A1 and A2 subfractions in toluene.²² The molecules are derived from those proposed by Acevedo et al.¹⁵ and were later modified to include oxygen.²²,²³ It was observed that the presence of hydrogen bonds is a key factor in aggregation and the structural connection between A1 and A2.²⁰,²²

In this paper, we investigate the behavior of asphaltene molecular models in solutions with different solvents: toluene (TOL), tetrahydrofuran (THF) and heptane (HEPT). For each solvent, we compute the size and the composition of the aggregates obtained after simulations at equilibrium from a starting point corresponding to dispersed asphaltene models and another starting point from a clustered configuration. Then, the driving forces of the aggregation mechanisms are discussed considering the relative weight of hydrogen bonds and dispersion interactions. Global properties of the obtained aggregates, such as the size and the dipole moments, are described. Finally, these results are compared to the size distribu-

tions observed by gel permeation chromatography separation followed by mass spectrometry analyses using an inductively coupled plasma for sample ionization: GPC-ICP-MS.¹⁶

# 2 Methodology

## 2.1 Structural models of asphaltenes

The structure of asphaltenes is one of the most controversial topics in petroleum chemistry. The controversy is centered on the organization of the aromatic cores within the asphaltene molecule (a single aromatic core or island against multiple aromatic cores or archipelago). The structures of the considered asphaltene molecular models were adopted to represent the A1 and A2 subfractions. The main structural characteristics of A1 and A2 subfractions consist of a lower H/C ratio and double bond equivalent (DBE) for A1. This makes molecules from the A1 subfraction more rigid than those from A2.²⁴ Acevedo et al.¹⁵ suggested two representative structures of asphaltene molecules which agree the above conditions. In this work, we considered the models used in a previous work,¹⁵,²² and two additional extreme cases were incorporated (see Figure 1). The first, A30 is a highly aromatic molecule derived from the A1-type of molecule, which increases the flatness characteristic and rigidity of the molecule. The second one, A40, is an archipelago-type molecule derived from an A2-type molecule, which has two aromatic nuclei connected by an aliphatic chain. These new models will help us to consider a small H/C ratio in the case of highly aromatic molecules and a high H/C ratio in the case of an archipelago-type molecule. Defining a set of molecular models to be representative of the asphaltene molecular family and used to implement molecular simulations may introduce a given bias in study like the current one. Heteroatoms, such as sulfur and aromatic nitrogen for example were added from pyridinic and thiophen structures. Although only a limited number of models were considered in this study and although they are
among the lightest of the observed asphaltene molecules,²⁵,²⁶ they were chosen to span the chemical diversity of the A1 and A2 subfractions keeping the implementation of the molecular simulations tractable.

The molecular models associated with the A1 and A2 subfractions are depicted in Figure 1. Table 1 summarizes the structural characteristics of the models and their corresponding labels along with the associated dipole moment obtained from an AM1-BCC calculation. Molecules A10 to A13 and A20 to A23 correspond to the structure considered in our previous work.²² Molecule A30 is a highly aromatic molecule, and molecule A40 is an archipelago-type molecule. Both models are proposed as extreme cases for subfractions A1 and A2, respectively. In order to implement the differences in the H/C ratio and the DBE between the A1 and A2 subfractions, the proposed A1-type molecules have two aliphatic cycles, while one cycle is exchanged with an n-octyl aliphatic chain in A2-type molecules. Figure S1 shows a general comparison of the molecules used in this work. This set of molecules will provide us with a way to explore the chemical space of the asphaltene fraction, including island and archipelago-type molecules, presence or absence of heteroatoms, various dipole moments, size of the aromatic nucleus and protic or aprotic groups. These molecules provide a way to tune the chemical interaction (hydrogen bonds and van der Waals interaction) and the rigidity of the molecule.

## 2.2 Simulation details

All simulations were carried out using the GROMOS96 force field with the 53A6 parameter set²⁷ using the GROMACS 2019.3 computer package.²⁸,²⁹ The main feature of this force field is its united-atom parameter (used for groups CH₃, CH₂ and CH₁), which allows an increase in the size of systems without sacrificing too much computational time. For each molecule, the force field parameters were assigned using parameters from the GROMOS force field, applying parameters from similar already-parametrized small organic molecules,

such as pyridine, thiophene, or phenol, which were used as building units. The charges were adapted using the AM1-BCC semiempirical method.³⁰,³¹ The Antechamber program within AmberTools²⁰³² was used to compute the partial charges of all molecules as suggested by Lemkul et al.³³ to produce GROMOS compatible force fields. The dipole moments obtained from the computed atomic charges are presented in Table 1. In the context of united atom force fields, atomic charges of CHₙ groups were summed up, as performed in reference 33. To describe hydrogen bonds, hydrogen atoms of protic groups such as phenol and carboxylic acids were considered explicitly.²⁰,²²

For each molecule, a simulation was performed in tetrahydrofuran, toluene, and heptane. Two types of simulation with different initial states were implemented, which are schematized in Figure 2. In the first case (A), the asphaltene molecules were randomly introduced into the box. In the second case (B), short molecular dynamics calculations were performed at a constant volume without solvent to promote the formation of an initial cluster including all the molecules. Then, this cluster obtained in vacuum was placed in a box filled by the considered solvents, and the asphaltene aggregates were obtained from its dissolution. This second scheme was implemented to provide a conceptual sketch of the experimental process associated with the dissolution of a powder containing the precipitated asphaltene aggregates into a given solvent before analysis.

The number of asphaltene molecules for each simulation was fixed at N=80 with an initial box size of 8×8×8 nm³ to obtain a concentration of ≈7wt.%. This number is within the range of the amount of molecules used in previous works.¹⁷,¹⁹,²⁰,²² In our previous work,²² we investigated the effect of the number of molecules used in the simulations and showed that 80 molecules is a converged value. In all simulations, periodic boundary conditions were implemented along with the particle mesh Ewalds method to manage electrostatic interactions.³⁴ A cutoff distance of 1.5 nm was considered for van der Waals interactions. All
<figure><img src="image_2.png" /><figcaption>Figure 1: Base molecular structure of the asphaltene molecules used in this study for the subfractions A1 and A2. In addition, the A3 and A4 molecules are presented.</figcaption></figure>
<table><caption>Table 1: Definition, labeling, and structural parameters of the molecules considered in this study. The base structure is illustrated in Figure 1. MW stands for molecular weight in gram per mole. The dipole moment was obtained from AM1-BCC calculation in order to adapt the partial charges of the GROMOS 54a7 force field.</caption><thead><tr><th></th><th>molecules</th><th>X<sub>1</sub></th><th>X<sub>2</sub></th><th>X<sub>3</sub></th><th>X<sub>4</sub></th><th>MW</th><th>DBE</th><th>H/C</th><th>μ (D)</th></tr></thead><tbody><tr><td rowspan="4">A1</td><td>A10</td><td>CH</td><td>CH</td><td>H</td><td>CH<sub>3</sub></td><td>572.836</td><td>23</td><td>1.000</td><td>1.98</td></tr><tr><td>A11</td><td>N</td><td>S</td><td>H</td><td>CH<sub>3</sub></td><td>591.870</td><td>23</td><td>0.976</td><td>5.09</td></tr><tr><td>A12</td><td>N</td><td>S</td><td>OH</td><td>CH<sub>3</sub></td><td>607.870</td><td>23</td><td>0.976</td><td>5.98</td></tr><tr><td>A13</td><td>N</td><td>S</td><td>H</td><td>COOH</td><td>621.854</td><td>24</td><td>0.928</td><td>3.95</td></tr><tr><td rowspan="4">A2</td><td>A20</td><td>CH</td><td>CH</td><td>H</td><td>CH<sub>3</sub></td><td>576.868</td><td>21</td><td>1.090</td><td>2.70</td></tr><tr><td>A21</td><td>N</td><td>S</td><td>H</td><td>CH<sub>3</sub></td><td>595.902</td><td>21</td><td>1.071</td><td>5.74</td></tr><tr><td>A22</td><td>N</td><td>S</td><td>OH</td><td>CH<sub>3</sub></td><td>611.902</td><td>21</td><td>1.071</td><td>6.71</td></tr><tr><td>A23</td><td>N</td><td>S</td><td>H</td><td>COOH</td><td>625.886</td><td>22</td><td>1.023</td><td>6.53</td></tr><tr><td>A3</td><td>A30</td><td>N</td><td>S</td><td>OH</td><td>CH<sub>3</sub></td><td>601.737</td><td>33</td><td>0.535</td><td>6.14</td></tr><tr><td>A4</td><td>A40</td><td>N</td><td>S</td><td>OH</td><td>CH<sub>3</sub></td><td>663.978</td><td>23</td><td>1.065</td><td>2.52</td></tr></tbody></table>

bonds were constrained using the LINCS algorithm.³⁵ The time step of the simulations was 2 fs, and the neighbor list update frequency was 10 fs. An initial minimization step was used to remove any high-energy structures using the conjugate gradient algorithm. The system was then allowed to reach thermal equilibrium at 300 K using the V-rescale thermostat³⁶ in the NVT ensemble for 200 ps. In a second equilibration step, the box density was adjusted in the NPT scheme with the V-rescale thermostat and the Parrinello-Rahman barostat³⁷,³⁸ at 1.0 bar pressure for 500 ps. Then, simulations in the NPT ensemble were performed at 1.0 bar and 300 K using the Parrinello-Rahman barostat and the Nose-Hoover thermostat.³⁹,⁴⁰ In scheme A, the initial configuration corresponds to a random dispersion of the asphaltene molecules, and the simulations were split into two steps. In accordance with the number of molecules included in the simulation, the equilibration time is long. A first step of 60 ns was performed to obtain aggregates and reach equilibrium. At this stage, dissociation and aggregation of molecules in large aggregates associated with an equilibrium state are observed in the simulations. Then, second step of 60 ns long was performed, which was considered as the production run and thus used in further

analyses to compute averages and perform the structural analyses of the aggregates. In the case of scheme B, starting from an aggregate state, only one step of 60 ns was implemented to analyze the stability of the aggregates in the different solvents.

## 2.3 Aggregate Analysis Methods

To perform the analyses of the aggregates obtained from the simulations, several structural parameters were considered and are illustrated in Figure 3. Molecules were considered to be in the same aggregate if the minimum distance between two atoms belonging to two different molecules (the distance of closest approach) was less than 3.5 Å. The analyses were performed using a homemade program that identifies the list and the composition of the molecular aggregates for each configuration extracted from the simulations. The aggregation number is the number of molecules that constitute an aggregate. This parameter is used as a guide to monitor both the equilibration and the expected average asphaltene aggregate size. The average aggregate number, $\overline{N}_{agg}$, was calculated from

$$ \overline{N}_{agg} = \frac{\sum_{i} i N_{i}}{\sum_{i} N_{i}} \qquad (1) $$
<figure><img src="image_3.png" /><figcaption>Figure 2: Scheme showing the construction of the initial configurations before molecular dynamics.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Structural parameters used in the simulations. a) Aggregation number is the number of molecules per aggregate, b) the hydrodynamic diameter, computed as the largest distance between two atoms in the same aggregate, c) the dipole moment of the overall aggregate.</figcaption></figure>

where $N_i$ is the number of aggregates containing $i$ molecules. The sums in Equation (1) extend from $i = 2$: that is, monomers were not included in the analysis.

Common ways to compute the aggregate size are either the largest distance between pairs of atoms within the same aggregate or the gyration radius.¹⁷,²² We showed in a previous work²² that the calculation of the aggregate size from the largest distance is more consistent with the hydrodynamic volume considered in gel permeation chromatography (GPC)²³,⁴¹ and was used as a size indicator in this work. This distance is recorded for each aggregate in each frame of

the simulation. These data are used to construct the histogram of aggregate sizes over all the simulations and provide information on the smallest and largest sizes of the aggregates. However, the size distribution obtained in this way is biased by the number of molecules considered in the simulations. Indeed, considering a total number of $N$ molecules, a total amount of $N/2$ dimers can be obtained, whereas only one aggregate of approximately $N$ molecules can be obtained. In order to take this bias into account, a distribution function $f(d)$ as a function of the hydrodynamic diameter $d$ of the aggregate can be constructed in the following way:

$$f(d) = \sum_{agg} \sum_{i \in n \ge 2}^{N} \frac{i\delta(d_i - d)}{N} \quad (2)$$

where $N$ is the total number of molecules in the system, $i$ is the number of molecules in an aggregate, and $d_i$ is the diameter of the current aggregate with $i$ molecules. The first sum runs over all aggregates obtained from all frames of the simulation, and the second sum runs over the number of molecules per aggregate, termed $i$ in the equation. Actually, the average size of aggregates containing $i$ molecules is weighted according to the number of aggre-
gates containing the same number of molecules that can be constructed by considering a total of N molecules: that is, $N/i$.

The dipole moment was computed for each aggregate throughout simulations using the atomic positions and the partial charges used in the force fields. This equation is additive and thus assumes that the polarizability of the solvent and the neighbor molecule is weak.

$$ \mu = \sum_{i=1}^{N} q_i r_i \qquad (3) $$

# 3 Results and discussion

## 3.1 Simulations in Heptane

As a first step, we considered the simulations in heptane. These results allow us to connect the aggregation behavior of the considered molecules with the behavior observed experimentally when considering the SARA classification and the A1 and A2 subfractions. Indeed, the asphaltene fraction is defined from the solubility of the molecule in heptane. Based on the definition of the aggregation number, when values close to the total number of molecules are reached, it indicates that the solubility of the molecules in the solvent is low, i.e., that the aggregate state is more stable than its dissociation.

<figure><img src="image_5.png" /><figcaption>Figure 4: Average aggregation number using heptane as solvent for all considered molecules. Circles on the left show an example of the two extreme cases obtained.</figcaption></figure>

Figure 4 compares the aggregation numbers of the molecules in heptane for all molecules considered in this work. The averages are computed

for the last 60 ns of the simulations following scheme A. The aggregation number commonly accepted in the Yen-Mullins paradigm is about at least 10 molecules per aggregates.⁸,⁹ Molecules A10, A11, A20 and A21 present an aggregation number lower than 10, which indicates that these molecules prefer to remain dispersed in the solvent. Molecules A10 and A20 are hydrocarbon molecules while molecules A11 and A21 contain few heteroatoms. However, any of them contain protic groups able to lead to the formation of hydrogen bonds. In that case, the only available driving force of the aggregation is the van der Waals interactions due to the aromatic nucleus or the dipole moment, which are rather low for these models; see Table 1. On the other hand, molecules A12, A13, A30, A22, A23 and A40 lead to an average aggregation number of 50 or more molecules. This situation is associated with the formation of large aggregates of at least 50 molecules on average. These molecules are thus most of the time included in large aggregates and can thus be considered as insoluble in heptane. These molecules include protic groups and are able to form hydrogen bonds that ease the formation of larger aggregates. In the following, we will consider only molecules that are insoluble in heptane, i.e. A12, A13, A22, A23, A30 and A40. This behavior agrees with the SARA classification of asphaltenes based on their solubility.

## 3.2 Simulations in THF and Toluene

This section compares the results of the simulations implemented in THF and toluene following the two schemes A and B, starting from an initially dispersed configuration or from a clustered configuration.

Figure 5 shows the histograms of the diameter of the aggregates obtained from the simulations of A12, A13, A30, A22, A23 and A40 molecules in toluene and THF following scheme B. The aggregate diameters are computed as the largest distance between two atoms belonging to the same aggregates. Only molecules with an average aggregation number larger than 20 were considered in this part. Indeed, as in-
Size distributions - Initially clustered

<figure><img src="image_6.png" /></figure>

Figure 5: Histograms of the diameter of the aggregates for molecules including protic groups (A12, A13, A30, A22, A23 and A40) starting from a clustered configuration. The histograms obtained for the same models from a dispersed configuration are available in Figure S2 in the supporting information. Aggregate size increases from right to left for further comparison with GPC experiments.

dicated in the previous section, these molecules are insoluble in heptane, which agrees with the definition of the asphaltene fraction in the SARA classification. The corresponding average aggregation numbers computed for these molecules in toluene and THF are presented in Figure S3 in the supporting information.

Regardless of the considered molecule, when comparing the aggregate sizes in toluene and THF, the largest aggregates are obtained in toluene. This result agrees with the values of the experimental solubility parameters of asphaltenes compared to those of toluene and THF. Indeed, the reported solubility parameters of asphaltene subfractions A1 and A2 are in the range of 20.5 to 22.7,⁴²,⁴³ which is closer to the solubility parameter of THF, 19.5, than to that of toluene 18.2.⁴⁴ The largest aggregates in toluene were thus expected and obtained from the simulations.

Considering scheme B, as shown in Figure 5 bottom panel, starting from a large aggregate, two different behaviors are obtained. In the case of molecules A12 and A22, the initial ag-

gregate is broken, and only small aggregates of 10 molecules or less, approximately 2.5 nm diameter, are formed in the simulation, as already obtained from scheme A. In contrast, from the simulations with molecules A13, A30, A23 and A40, the aggregate size distribution displays a bimodal shape associated with the presence of aggregates with a diameter larger than 7 nm and on average more than 20 molecules per aggregate (see Figure S3 in the supporting information). These results suggest that in the case of molecules with structural properties such as A13, A30, A23 and A40, including heteroatoms and free carboxylic groups, when initially packed aggregates are placed in solution, THF and toluene solvents are not strong enough to disrupt the aggregates and dissolve molecules or small aggregates. This information is of highest interest considering the experimental setup for asphaltene aggregate analyses in GPC-ICP-MS experiments where asphaltene fractions are first precipitated and then analyzed after dissolution in THF.²³,⁴⁵⁻⁴⁷ This indicates that under certain conditions, the ag-
gregates analyzed experimentally by GPC-ICP-MS already existed in the crude and were preserved in THF.²³,⁴⁵⁻⁴⁷ The comparison with experimental results will be discussed in the last section of the paper.

In the case of the initially dispersed systems (see Figure S2 in the supporting information), systems A12, A13, A22 and A23 exhibit a similar behavior. Only a few molecules aggregate, and the average aggregation numbers are in the range between 5 to at most 10 molecules per aggregate, with the smallest values obtained for THF. Again, these results in toluene agree with the definition of the asphaltene fraction according to its solubility in the SARA classification, assuming that aggregates are too small to precipitate and are thus soluble in toluene. These molecules are thus slightly more soluble in THF than in toluene. These results support our hypothesis that aggregates observed experimentally are directly those that already exist in the crude and not formed during the analysis.

Going further in the analyses of Figure 5, the results may be linked to the chemical characteristics of the models and the inter-molecular energies displayed figure S10. Molecules A12 and A22 contain nitrogen and sulfur heteroatoms along with an OH group attached to the nucleus. Although their dipole moment is strong, the formation of large aggregates from the formation of hydrogen bonds including the OH groups is hindered due to steric effects associated with the buried position of this protic group. Moreover, this OH group bound easily with THF molecules, which are small enough to form hydrogen bonds. One can see on figure S10, that this molecules are those which make the largest number of hydrogen bonds with THF molecules. These phenomena limit the aggregate size both in toluene and in THF.

Molecules A13 and A23 contain a carboxyl group bonded to a terminal carbon atom of an alkyl chain, and large aggregates were formed during the simulations in THF. This is in agreement with the results from other work²⁰ indicating that the presence of heteroatoms in the terminal positions of side chains enhances a more effective self-assembly compared to the aromatic ring positions.⁴⁸,⁴⁹ This ability is di-

rectly related to the formation of hydrogen bonds between asphaltene molecules. One can see on figure S10, that this molecules are those which make the largest number of hydrogen bonds between asphaltene molecules. But also, a large number of hydrogen bonds with THF molecules is obtained. In consequence, the size of aggregates is significantly smaller when following scheme A starting from a dispersed situation. When the system is initially dispersed, THF molecules act as competitive inhibitors and limit the ability of asphaltene molecules to form hydrogen bonds that promote π-stacking-type interactions. This leads to the formation of smaller aggregates.

In the case of the A30 molecule (a highly aromatic molecule) and A40 molecule (archipelago type), whatever the starting configuration, and following scheme A or B, large aggregates with more than 60 molecules are obtained with a diameter larger than 10 nm. This reveals the low solubility of these molecules in toluene or THF. These models were considered to include in the discussion a flat and rigid molecule (A30) and a highly flexible molecule (A40). Representative aggregates obtained from simulations are depicted in Figure 6. In the case of the A30 molecule (rigid, flat and with a large aromatic nucleus) the π-stacking interaction between molecules is clearly favored. This leads to the formation of tubular clusters that can be aggregated in large structures by connecting the cylinders due to side chain interactions, such as hydrogen bonds, or T-shaped interactions between two aromatic nuclei. From figure S10 again, A30 molecules form only a few amount of hydrogen bonds indicating that π-stacking interactions are favored. The combination of π-stacking interactions and side chain interactions leads to the formation of the largest aggregates. Representative structures of aggregates obtained from simulations of A30 are depicted in Figure S4. Figures S5 to S8 present the computed values of the shape anisotropy of the aggregates, which provides an idea of the overall shape of the aggregates.¹⁷,²²

In the case of molecule A40, of the archipelago type, large aggregates are obtained from the agglomeration of small globular aggregates, as dis-
<figure><img src="image_7.png" /><figcaption>Figure 6: Representative structures of aggregates obtained from simulations with molecule A30 on the left and molecule A40 on the right illustrating the various types of interactions.</figcaption></figure>

played in Figure 6. These aggregates combine small π-stacking interactions or T-shape interactions with strong hydrogen bond interactions between side chains allowed by the high flexibility of the structure. This leads to the formation of large aggregates. In a recent work, Yang et al.⁵⁰ suggested a similar conclusion and propose that when asphaltenes become more unstable, an increase in the amount of π-stacking interactions is observed, but this is concomitant with the formation of disordered structure of aggregates due to side chain interactions.

## 3.3 Link between simulations and chromatography results

GPC-ICP-MS analyses of the asphaltene fraction of crude oil are a common technique to characterize such complex matrices.¹⁶,²³ The separation of the compounds of the sample is based on their size, with the larger compounds requiring a longer time to pass through the column.²³,⁴⁵ The presence of asphaltene aggregates is detected from the presence of vanadyl porphyrins in the crude oil which are embedded or linked to the aggregates.²³,⁵¹,⁵² In this section, we show that it is possible to make a conceptual comparison between the GPC-ICP-MS profile and the size distribution of the ag-

gregates obtained from the molecular dynamics simulations. We chose molecules A12 and A22 to illustrate this comparison because they belong to the A1 and A2 families and their solubility behaviors are representative of these fractions: they are insoluble in heptane and soluble in toluene. This choice is a strong limitation, but the aim here is to illustrate the capability of the coupling between the GPC techniques and the molecular dynamics methodology to work in synergy.

Figure 7 b) shows the chromatograms obtained after GPC-ICP-MS analyses performed by our group in a previous work²³ on a sample of Hamaca crude oil. A complete description of the experimental setup is provided in the work of Putman et al.⁴¹. Briefly, in order to perform the separation, three polymeric Shodex preparative GPC columns were connected in series with a pore size of 50, 300 and 1500 Å. Three samples were analyzed: the complete asphaltene fraction and then the A1 and A2 subfractions.²³,⁴⁶,⁴⁷ In this sample, it was determined that the relative proportions of A1 and A2 were 40% and 60%, respectively. The graph shows the variation in the signal intensity of the detection of the ⁵¹V isotope as a function of elution time, which is a conventional way to follow asphaltene aggregates. The intensity of the sig-
<figure><img src="image_8.png" /><figcaption>Figure 7: a) Size distributions in THF from simulations of molecules A12 and A22 following scheme B with initially clustered systems. The combination of A12 and A22 distributions was performed using the experimental proportions of A1 and A2 of Hamaca crude oil, that is, $A1 \times 0.4 + A2 \times 0.6$.²³ b) GPC-ICP-MS chromatograms measured for the asphaltene fractions A1 and A2 of a Hamaca crude oil.²³ Super high molecular weight (SHMW) region, high molecular weight (HMW) region, medium molecular weight (MMW) and low molecular weight (LMW) regions are labeled and delimited by vertical dashed lines.</figcaption></figure>

nal is proportional to the number of molecules, while the elution time is associated with the size or molar mass of the molecules. We can observe the presence of four regions, labeled super high molecular weight (SHMW), high molecular weight (HMW), medium molecular weight (MMW) and low molecular weight (LMW). The complete experimental setup is detailed in a previous work 23. After precipitation of the asphaltene fraction and A1 and A2 subfractions, the powder containing aggregates was dissolved in THF and injected into the apparatus using THF as the eluting solvent.

In this case, we followed scheme B, considering the simulations starting from a clustered initial configuration further dissolved in THF. Indeed, this scheme can be associated with the dissolution of the precipitated asphaltene in THF before injection into the apparatus to provide us with a conceptual sketch of the experimental process. The size distributions of the aggregates obtained from simulations considering A12 and A22 in THF and following scheme B were computed using equation 2. They are depicted in figure 7 a) along with a reconstructed GPC-ICP-MS profile using the relative proportions of A1 and A2 in the experimental sample. First, it is observed that the shape distributions are similar, differing in that the low molecular weight (LMW) aggregates of A22 are larger than those of A12 by 1 nm, and A22 shows a greater intensity in the high molecular weight (HMW) region. From this comparison, it is possible to generate a scheme of the aggregate formation during the experimental GPC-ICP-MS measurement of subfractions A1 and A2.

Although OH groups are present in molecules A12 and A22, once dissolved in THF, the initial cluster is split into smaller and more stable aggregates. According to this result, the clearest interpretation is that the initial large aggregates are split into smaller fractions, and these fractions remain in equilibrium or are metastable in the system. This interpretation is consistent with Cheng's model,⁵³ where asphaltene undergoes fractionation into nanoaggregates of asphaltenes, which can be interpreted as a top-down aggregate formation different from the down-to-top association proposed in the micel-
lar model of Gray et al.¹⁰. A prominent feature of the aggregate structure is the acquired shape. Representative aggregates from simulations of A12 and A22 are depicted in Figure 7 a). Whatever the size, the π-stacking interaction prevails in the aggregates of A12 molecules, leading to a dense structure. In contrast, in the aggregates of A22 molecules, a less dense structure is observed due to the aliphatic chain of the A2 family, and more hydrogen bonds are involved in building large aggregates.

Additional information can be obtained from the analyses of molecular dynamics simulations by considering the number of molecules by aggregate in each weight region. Figure 8 shows a box-plot representation of the distribution of the number of molecules per aggregate in each molecule weight region of the chromatogram. It is interesting to note that by using the results of the simulations, we can roughly estimate both the size of the aggregates and the number of molecules per aggregate in each region. On average, SHMW aggregates are composed of 60 molecules, while the value decreases to 20 for HMW aggregates, and finally, fewer than 5 molecules are enough to obtain MMW aggregates. The LMW region is limited by the choice of asphaltene molecule. Indeed, the mass values corresponding to this region range from 162 to 1,100 Da²³ corresponding to the monomers in our case.

<figure><img src="image_9.png" /><figcaption>Figure 8: Number of molecules present per molecular weight region observed for these models.</figcaption></figure>

## 3.4 Dipole moment analyses

In this last section, we finally analyzed the link between the dipole moment of the aggregates and their size or the mass region to which they belong. Figure 9 shows the dipole moment of the aggregates obtained during simulation as a function of time and as a function of the number of molecules per aggregate. The colors represent the weight regions of the aggregates. Examining the dipole moment values of the individual molecules reported in table 1, the dipole moment of aggregates is on average lower than the value of the individual molecule, except for the A40 molecule. In a few cases, such as A13, A30 or A40, a discontinuity is obtained because no aggregates were observed throughout the simulation with those numbers of molecules. Whatever the considered molecule, there is a clear trend indicating that increasing the number of molecules per aggregate, and thus the size of the aggregate, leads to a decrease in the dipole moment of the whole aggregate. This behavior is typical of systems grown randomly, leading to the partial cancellation of the individual dipole moments of the molecules composing the aggregates. These results enforce the conclusion that to obtain large aggregates, the formation of π-stacking interactions is concomitant with side chain interactions, thus leading to more disordered structures.

In several cases, the values of the dipole moment present a wide dispersion for a given aggregate size. Whatever the molecule, this is true for small aggregates for which the combination of a small number of molecules can lead to an addition or a subtraction of their individual dipole moment. In particular, for molecule A40, the high flexibility of this molecule leads to a variance of the dipole moment in the LMW region, which is very high. However, as with other molecules, increasing the size of the aggregates ultimately leads to a decrease in the whole dipole moment of the aggregates. Indeed, aggregates obtained from flexible molecules adopt a more globular structure for which, during the simulation, the size, shape and relative arrangement of molecules present large fluctuations during the simulations. Figures S5, S6,
<figure><img src="image_10.png" /><figcaption>Figure 9: Dipole moments of aggregates obtained for molecules A12, A13, A30, A22, A23 and A40 grouped by weight regions. From left to right, the plots represent the dipole moment as a function of time, as a function of the molecules per aggregate, and finally as a histogram for the simulations following scheme B.</figcaption></figure>
S7 and S8 in the supporting information illustrate that point considering the dispersion of the computed shape anisotropy.

In the case of highly aromatic systems, previous works indicated that high dipole moment values can be attributed to a small distance between stacking layers. This small distance is a consequence of a larger aromatic ring structure and a larger average number of aligned stacking layers.⁵⁴ In the case of A30 molecules, the rigidity of the molecule leads to the formation of these rod-like structures that produce an intense dipole moment. However, large aggregates are obtained by the random cohesion of these rod-like structures, leading to the cancellation of the dipole moment (see Figure 6).

## 4 Conclusions

In this paper, considering a set of asphaltene molecular models representative of A1 and A2 subfractions of asphaltene, we investigated the aggregation of asphaltene molecules in heptane, toluene and THF solvents using molecular dynamics simulations. The asphaltene molecules were chosen to tune the structural characteristics of the subfractions (H/C ratio, DBE) along with the possible intermolecular interactions and the rigidity of the molecules. Two schemes were followed, starting either from a dispersed configuration of the asphaltene molecules or from a clustered configuration. The asphaltene molecules show an acceptable response to the expected experimental behavior. The solubility in toluene and heptane is consistent with the definition of asphaltenes from SARA classification. The results show that the formation of large aggregates requires the concomitant engagement of π-stacking and hydrogen bonding interactions. Implementing molecular dynamics simulations starting from an initial clustered configuration of the molecules provides a conceptual sketch of the experimental process of asphaltene sample analysis. We showed that THF is able to preserve large aggregates after dissolution, which validates the hypothesis that experimentally analyzed aggregates are not formed in a bottom-up sketch during the experiment-

tal setup but are obtained from the fragmentation of larger aggregates initially present in the crude. Finally, we link the experimental GPC-ICP-MS profile and the size distribution of the aggregates obtained from the simulations. This allows us to characterize asphaltene aggregates belonging to the different weight regions in terms of size, structure, number of molecules and dipole moment. We showed that the overall dipole moment decreases in accordance with the number of molecules due to randomness or the globular structure of the large aggregates.

**Acknowledgement** This work was achieved using HPC resources from GENCI-CINES (grant 2021-A0100806920), MCIA (Mésocentre de Calcul Intensif Aquitain) and the Université de Pau et des Pays de l'Adour.

## Supporting Information Available

Supporting information is available and provides (i) complementary information about molecular models and (ii) additional data on asphaltene aggregates, including all models for completeness and (iii) the force field parameters (itp files) of the asphaltene molecules and the solvent molecule.

## References

(1) Fan, T.; Buckley, J. S. Rapid and Accurate SARA Analysis of Medium Gravity Crude Oils. *Energy & Fuels* **2002**, *16*, 1571–1575.

(2) Haji-Akbari, N.; Masirisuk, P.; Hoepfner, M. P.; Fogler, H. S. A Unified Model for Aggregation of Asphaltenes. *Energy and Fuels* **2013**, *27*, 2497–2505.

(3) Chacón-Patiño, M. L.; Rowland, S. M.; Rodgers, R. P. Advances in Asphaltene Petroleomics. Part 1: Asphaltenes Are Composed of Abundant Island and Archipelago Structural Motifs. *Energy & Fuels* **2017**, *31*, 13509–13518.
(4) Dutta Majumdar, R.; Montina, T.; Mullins, O. C.; Gerken, M.; Hazendonk, P. Insights into asphaltene aggregate structure using ultrafast MAS solid-state 1H NMR spectroscopy. *Fuel* **2017**, *193*, 359–368.

(5) Leontaritis, K. J.; Kawanaka, S.; Mansoori, G. A. Descriptive accounts of thermodynamic and colloidal models of asphaltene flocculation. *Proceedings of Symposium on Thermodynamic and Transport Properties of Fluids*. **1988**, *33*, 196–204.

(6) Mullins, O. C.; Betancourt, S. S.; Cribbs, M. E.; Dubost, F. X.; Creek, J. L.; Andrews, A. B.; Venkataramanan, L. The Colloidal Structure of Crude Oil and the Structure of Oil Reservoirs. *Energy & Fuels* **2007**, *21*, 2785–2794.

(7) Headen, T. F.; Boek, E. S.; Skipper, N. T. Evidence for Asphaltene Nanoaggregation in Toluene and Heptane from Molecular Dynamics Simulations. *Energy & Fuels* **2009**, *23*, 1220–1229.

(8) Mullins, O. C. The Modified Yen Model. *Energy & Fuels* **2010**, *24*, 2179–2207.

(9) Mullins, O. C.; Sabbah, H.; Eyssautier, J.; Pomerantz, A. E.; Barré, L.; Andrews, A. B.; Ruiz-Morales, Y.; Mostowfi, F.; McFarlane, R.; Goual, L.; Lepkowicz, R.; Cooper, T.; Orbulescu, J.; Leblanc, R. M.; Edwards, J.; Zare, R. N. Advances in Asphaltene Science and the Yen–Mullins Model. *Energy & Fuels* **2012**, *26*, 3986–4003.

(10) Gray, M. R.; Tykwinski, R. R.; Stryker, J. M.; Tan, X. Supramolecular Assembly Model for Aggregation of Petroleum Asphaltenes. *Energy & Fuels* **2011**, *25*, 3125–3134.

(11) Akbarzadeh, K.; Bressler, D. C.; Wang, J.; Gawrys, K. L.; Gray, M. R.; Peter K. Kilpatrick,; Yarranton§, H. W. Association Behavior of Pyrene Compounds as Models for Asphaltenes†. *Energy & Fuels* **2005**, *19*, 1268–1271.

(12) Alshareef, A. H. Asphaltenes: Definition, Properties, and Reactions of Model Compounds. *Energy & Fuels* **2019**, *34*, 16–30.

(13) Acevedo, S.; Castro, A.; Vásquez, E.; Marcano, F.; Ranaudo, M. A. Investigation of Physical Chemistry Properties of Asphaltenes Using Solubility Parameters of Asphaltenes and Their Fractions A1 and A2. *Energy & Fuels* **2010**, *24*, 5921–5933.

(14) Gutiérrez, L. B.; Ranaudo, M. A.; Bernardo Méndez,; Acevedo, S. Fractionation of Asphaltene by Complex Formation with p-Nitrophenol. A Method for Structural Studies and Stability of Asphaltene Colloids. *Energy & Fuels* **2001**, *15*, 624–628.

(15) Acevedo, S.; Castillo, J.; Vargas, V.; Castro, A.; Delgado, O.; Cortés, F. B.; Franco, C. A.; Bouyssiere, B. Suppression of Phase Separation as a Hypothesis to Account for Nuclei or Nanoaggregate Formation by Asphaltenes in Toluene. *Energy & Fuels* **2018**, *32*, 6669–6677.

(16) Castillo, J.; Gonzalez, G.; Bouyssiere, B.; Vargas, V. Asphaltenes, subfractions A1 and A2 aggregation and adsorption onto RH-SiO2 nanoparticles: Solvent effect on the aggregate size. *Fuel* **2023**, *331*, 125635.

(17) Headen, T. F.; Boek, E. S.; Jackson, G.; Totton, T. S.; Müller, E. A. Simulation of Asphaltene Aggregation through Molecular Dynamics: Insights and Limitations. *Energy & Fuels* **2017**, *31*, 1108–1125.

(18) Law, J. C.; Headen, T. F.; Jiménez-Serratos, G.; Boek, E. S.; Murgich, J.; Müller, E. A. Catalogue of Plausible Molecular Models for the Molecular Dynamics of Asphaltenes and Resins Obtained from Quantitative Molecular Representation. *Energy & Fuels* **2019**, *33*, 9779–9795.
(19) Silva, H. S.; Sodero, A. C. R.; Bouyssiere, B.; Carrier, H.; Korb, J.-P.; Alfarra, A.; Vallverdu, G.; Bégué, D.; Baraille, I. Molecular Dynamics Study of Nanoaggregation in Asphaltene Mixtures: Effects of the N, O, and S Heteroatoms. *Energy & Fuels* **2016**, *30*, 5656–5664.

(20) Santos Silva, H.; Alfarra, A.; Vallverdu, G.; Bégué, D.; Bouyssiere, B.; Baraille, I. Impact of H-Bonds and Porphyrins on Asphaltene Aggregation As Revealed by Molecular Dynamics Simulations. *Energy & Fuels* **2018**, *32*, 11153–11164.

(21) Santos Silva, H.; Alfarra, A.; Vallverdu, G.; Bégué, D.; Bouyssiere, B.; Baraille, I. Asphaltene aggregation studied by molecular dynamics simulations: role of the molecular architecture and solvents on the supramolecular or colloidal behavior. *Petroleum Science* **2019**, *16*, 669–684.

(22) Villegas, O.; Salvato Vallverdu, G.; Bouyssiere, B.; Acevedo, S.; Castillo, J.; Baraille, I. Molecular Cartography of A1 and A2 Asphaltene Subfractions from Classical Molecular Dynamics Simulations. *Energy & Fuels* **2020**, *34*, 13954–13965.

(23) González, G.; Acevedo, S.; Castillo, J.; Villegas, O.; Ranaudo, M. A.; Guzmán, K.; Orea, M.; Bouyssiere, B. Study of Very High Molecular Weight Cluster Presence in THF Solution of Asphaltenes and Subfractions A1 and A2, by Gel Permeation Chromatography with Inductively Coupled Plasma Mass Spectrometry. *Energy & Fuels* **2020**, *34*, 12535–12544.

(24) Pekerar, S.; Lehmann, T.; Méndez, B.; Acevedo, S. Mobility of Asphaltene Samples Studied by 13C NMR Spectroscopy. *Energy & Fuels* **1999**, *13*, 305–308.

(25) Schuler, B.; Meyer, G.; Peña, D.; Mullins, O. C.; Gross, L. Unraveling the

Molecular Structures of Asphaltenes by Atomic Force Microscopy. *Journal of the American Chemical Society* **2015**, *137*, 9870–9876.

(26) Schuler, B.; Fatayer, S.; Meyer, G.; Rogel, E.; Moir, M.; Zhang, Y.; Harper, M. R.; Pomerantz, A. E.; Bake, K. D.; Witt, M.; Peña, D.; Kushnerick, J. D.; Mullins, O. C.; Ovalles, C.; van den Berg, F. G. A.; Gross, L. Heavy Oil Based Mixtures of Different Origins and Treatments Studied by Atomic Force Microscopy. *Energy & Fuels* **2017**, *31*, 6856–6861.

(27) Oostenbrink, C.; Villa, A.; Mark, A. E.; Van Gunsteren, W. F. A biomolecular force field based on the free enthalpy of hydration and solvation: The GROMOS force-field parameter sets 53A5 and 53A6. *Journal of Computational Chemistry* **2004**, *25*, 1656–1676.

(28) Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, flexible, and free. *Journal of Computational Chemistry* **2005**, *26*, 1701–1718.

(29) Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; Hess, B.; Lindahl, E. GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit. *Bioinformatics* (Oxford, England) **2013**, *29*, 845–54.

(30) Jakalian, A.; Bush, B. L.; Jack, D. B.; Bayly, C. I. Fast, Efficient Generation of High-Quality Atomic Charges. AM1-BCC Model: I. Method. *Journal of Computational Chemistry* **2000**, *21*, 132–146.

(31) Jakalian, A.; Jack, D. B.; Bayly, C. I. Fast, efficient generation of high-quality atomic charges. AM1-BCC model: II. Parameterization and validation. *Journal of Computational Chemistry* **2002**, *23*, 1623–1641.
(32) Case, D. A.; Cheatham, T. E.; Darden, T.; Gohlke, H.; Luo, R.; Merz, K. M.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J. The Amber biomolecular simulation programs. *Journal of Computational Chemistry* **2005**, *26*, 1668–1688.

(33) Lemkul, J. A.; Allen, W. J.; Benan, D. R. Practical Considerations for Building GROMOS-Compatible Small-Molecule Topologies. *Journal of Chemical Information and Modeling* **2010**, *50*, 2221–2235.

(34) Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: An N log(N) method for Ewald sums in large systems. *The Journal of Chemical Physics* **1993**, *98*, 10089–10092.

(35) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A linear constraint solver for molecular simulations. *Journal of Computational Chemistry* **1997**, *18*, 1463–1472.

(36) Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. *The Journal of Chemical Physics* **2007**, *126*, 014101.

(37) Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. *Journal of Applied Physics* **1981**, *52*, 7182–7190.

(38) Nosé, S.; Klein, M. Constant pressure molecular dynamics for molecular systems. *Molecular Physics* **1983**, *50*, 1055–1076.

(39) Nosé, S. A molecular dynamics method for simulations in the canonical ensemble. *Molecular Physics* **1984**, *52*, 255–268.

(40) Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions. *Physical Review A* **1985**, *31*, 1695–1697.

(41) Putman, J. C.; Gutiérrez Sama, S.; Barrère-Mangote, C.; Rodgers, R. P.; Lobinski, R.; Marshall, A. G.;

Bouyssière, B.; Giusti, P. Analysis of Petroleum Products by Gel Permeation Chromatography Coupled Online with Inductively Coupled Plasma Mass Spectrometry and Offline with Fourier Transform Ion Cyclotron Resonance Mass Spectrometry. *Energy & Fuels* **2018**, *32*, 12198–12204.

(42) Hansen, C. M. *Hansen Solubility Parameters: A User's Handbook*, 2nd ed.; CRC press, 2007.

(43) Acevedo, S.; Castro, A.; Vásquez, E.; Marcano, F.; Ranaudo, M. A. Investigation of Physical Chemistry Properties of Asphaltenes Using Solubility Parameters of Asphaltenes and Their Fractions A1 and A2. *Energy & Fuels* **2010**, *24*, 5921–5933.

(44) Hildebrand, J. H.; Scott, R. L. *The solubility of nonelectrolytes*; Reinhold Pub. Corp.: New York, 1950; OCLC: 1137083.

(45) Gascon, G.; Vargas, V.; Feo, L.; Castellano, O.; Castillo, J.; Giusti, P.; Acevedo, S.; Lienemann, C. P.; Bouyssiere, B. Size Distributions of Sulfur, Vanadium, and Nickel Compounds in Crude Oils, Residues, and Their Saturate, Aromatic, Resin, and Asphaltene Fractions Determined by Gel Permeation Chromatography Inductively Coupled Plasma High-Resolution Mass Spectrometry. *Energy & Fuels* **2017**, *31*, 7783–7788.

(46) Gascon, G.; Negrín, J.; García-Montoto, V.; Acevedo, S.; Lienemann, C. P.; Bouyssiere, B. Simplification of Heavy Matrices by Liquid-Liquid Extraction: Part i - How to Separate LMW, MMW, and HMW Compounds in Maltene Fractions of V, Ni, and S Compounds. *Energy & Fuels* **2019**, *33*, 1922–1927.

(47) Gascon, G.; Negrín, J.; Montoto, V. G.; Acevedo, S.; Lienemann, C. P.; Bouyssiere, B. Simplification of Heavy
Matrices by Liquid-Solid Extraction: Part II-How to Separate the LMW, MMW, and HMW Compounds in Asphaltene Fractions for V, Ni, and S Compounds. *Energy & Fuels* **2019**, *33*, 8110–8117.

(48) Zhang, L. L.; Yang, G. H.; Wang, J. Q.; Li, Y.; Li, L.; Yang, C. H. Study on the polarity, solubility, and stacking characteristics of asphaltenes. *Fuel* **2014**, *128*, 366–372.

(49) Ekramipooya, A.; Mirzaee Valadi, F.; Farisabadi, A.; Gholami, M. R.; Valadi, F. M.; Farisabadi, A.; Gholami, M. R. Effect of the heteroatom presence in different positions of the model asphaltene structure on the self-aggregation: MD and DFT study. *Journal of Molecular Liquids* **2021**, *334*, 116109.

(50) Yang, Y.; Headen, T. F.; Zhang, J.; Hoepfner, M. P. Solvent Effects on the Structure of Petroleum Asphaltenes. *Energy & Fuels* **2021**, *35*, 13743–13755.

(51) Putman, J. C.; Moulian, R.; Barrère-Mangote, C.; Rodgers, R. P.; Bouyssiere, B.; Giusti, P.; Marshall, A. G. Probing Aggregation Tendencies in Asphaltenes by Gel Permeation Chromatography. Part 1: Online Inductively Coupled Plasma Mass Spectrometry and Offline Fourier Transform Ion Cyclotron Resonance Mass Spectrometry. *Energy & Fuels* **2020**, *34*, 8308–8315.

(52) Chacón-Patiño, M. L.; Moulian, R.; Barrère-Mangote, C.; Putman, J. C.; Weisbrod, C. R.; Blakney, G. T.; Bouyssiere, B.; Rodgers, R. P.; Giusti, P. Compositional Trends for Total Vana-dium Content and Vanadyl Porphyrins in Gel Permeation Chromatography Fractions Reveal Correlations between Asphal-tene Aggregation and Ion Production Effi-ciency in Atmospheric Pressure Photoion-ization. *Energy & Fuels* **2020**, *34*, 16158–16172.

(53) Cheng, B.; Zhao, J.; Yang, C.; Tian, Y.; Liao, Z. Geochemical Evolution of Oc-cluded Hydrocarbons inside Geomacro-molecules: A Review. *Energy & Fuels* **2017**, *31*, 8823–8832.

(54) Farooq, U.; Laedre, S.; Gawel, K. Review of Asphaltenes in an Electric Field. *Energy & Fuels* **2021**, *35*, 7285–7304.