Influence of Temperature on the Adsorption and Diffusion of Heavy Oil in Quartz Nanopore: A Molecular Dynamics Study

Abstract

The desorption of heavy oil is one of the important indicators affecting the development efficiency of the remaining oil in nanopores. However, the study of the adsorption and diffusion mechanisms of heavy oil molecules in nanopores remains scarce. In this work, the influences of temperature on the adsorption and diffusion properties of the heavy oil four-fractions in quartz nanopore have been investigated via molecular dynamics simulations. Our results show that the heavy oil molecules will form a denser multilayer adsorption oil layer on the nanopore surface, and high temperature can alter the adsorption behaviors of the heavy oil four-fractions. As the temperature increases, the saturate molecules are desorbed from the nanopore surfaces, but the aromatic, resin, and asphaltene molecules maintain a tendency to aggregate towards the nanopore surface. In particular, the agglomeration behaviors of most saturate, aromatic and asphaltene molecules in nanopore can be suppressed by the confined space compared with the heavy oil molecules in oil droplet. In addition, the influence of temperature on the movement of heavy oil molecules in nanopore decreases compared with the oil molecules in a heavy oil droplet due to the confined space and adsorption effect. Interestingly, there is a competition phenomenon between the adsorption and diffusion of aromatic, resin, and asphaltene molecules in the nanopore, resulting in different adsorption behaviors with the increase in temperature. The results obtained in this paper will provide molecular-level theoretical guidance for understanding the adsorption and desorption mechanisms of heavy oil in nanopores.

1. Introduction

The continuous consumption of conventional oil resources is leading to the intensification of the energy crisis, thus unconventional oil reservoir resources such as heavy oil, extra-heavy oil, and bitumen will play a crucial role in the global energy supply. Studies show that heavy oil and bitumen are predicted to make up 70% of the total remaining oil reservoirs [1]. However, the physical properties of heavy oil, such as high viscosity, high density, and poor fluidity, seriously hinder the recovery efficiency [2]. Nowadays, thermal recovery is one of the main methods of heavy oil recovery, which can effectively reduce the viscosity of heavy oil by heating it. Nevertheless, the growth rate of heavy oil recovery will slow down and gradually stabilized during the middle and later stages of thermal recovery [3] due to the significant reduction of thermal efficiency and the adsorption effect of heavy oil on rock surfaces. Over the years, there have been numerous efforts to improve the techniques of enhanced oil recovery (EOR). For instance, Li et al. [4] found that a hot solvent with a higher carbon number can result in a larger latent heat of vaporization and a better viscosity reduction effect by dissolution and heating. Their results show that the solvent vapor can render greater sweep efficiency (75%), displacement efficiency (85%), and high oil recovery (65%). In addition, Wang et al. [5] found that the flue gas can promote the steam chamber to develop upward and then expand laterally when steam and flue gas are mixed injection, which greatly improves the steam chamber sweep coefficient and the recovery degree. Similarly, Fan et al. [6] carried out high temperature and pressure experiments and different huff and puff experiments of the sand pack. It was found that the multi-gas assisted steam huff and puff process can increase the recovery by 2–5%. According to the above research, although extensive investigations have been carried out on the EOR techniques, there is still a lack of deep understanding of the mechanisms for enhancing remaining heavy oil recovery.

The adsorption phenomena of heavy oil in rock cracks [7] and pipelines [8] are ubiquitous, which seriously hinders the deep development and utilization of heavy oil resources. Therefore, the occurrence and adsorption states of heavy oil on rock surfaces are vital factors to further improve the EOR techniques, which have drawn growing research interests recently. For instance, Huang et al. [9] probed the interaction behaviors of two distinct petroleum products, diluted bitumen and conventional crude oil, with molecularly smooth mica surfaces in toluene solutions using a surface forces apparatus and an atomic force microscope. Their results indicate that the adsorption and interaction between two types of oil on the mica surface are highly dependent on concentration and adsorption time. In addition, Mohammadi et al. [10] investigated the adsorption behaviors of six asphaltene samples of various origins onto the magnetite surface by experimental measures and adsorption isotherm models, which indicate that an increase in the nitrogen content and aromatic nature of asphaltenes increased their adsorption on magnetite. As mentioned above, although the adsorption strength of heavy oil can be characterized by the sophisticated apparatus, it is still difficult to directly observe the underlying adsorption mechanisms of heavy oil due to the limited technologies and devices.

Molecular dynamics (MD) simulation is one of the important means to investigate the molecular-level physical mechanisms, which has been extensively used to study the adsorption behaviors of organic matter. Particularly, the investigation of adsorption mechanisms of heavy oil has received extensive attention in recent years [8,11,12,13,14]. For instance, Lyu et al. [8] studied the effect of the heavy oil composition and molecular structures on adhesion via experiments and MD simulations. Their results show that the addition of asphaltene increases the adhesion of heavy oil by as much as 16.60–83.35%. In addition, Bablu Meghwal et al. [11] found a dense interfacial layer adjacent to the rock surface by MD simulations, which in turn governs the macroscopic adsorption behavior between the heavy oil and rock. Moreover, Ji et al. [12] studied the adsorption of heavy oil droplets on different hydrophilic and hydrophobic silica surfaces by MD simulations, and found that the resin and asphaltene molecules on different surfaces showed different aggregated structures. Particularly, the oil droplet is easier to adsorb on the fully hydrophobic system due to the face-to-face stacking interaction of resin and asphaltene molecules. Zhang et al. [13] investigated the adsorption and distribution of heavy oil subfractions on an Na-montmorillonite surface via a similar simulation method, which can observe the adsorption zone and free zone of heavy oil. Liu et al. [14] studied the effects of molecular polarity on the adsorption and desorption behaviors of asphaltenes on silica surfaces via the experiments and MD simulations, and their results show that the asphaltenes exhibit different adsorption behaviors on surfaces with different wettability. In addition, the influence of solvent [15,16] and surfactant [17] on the adsorption behaviors of heavy oil has also been studied by some groups via MD simulations. However, despite the importance of adsorption behaviors of heavy oil being widely researched, the adsorption mechanisms of heavy oil four-fractions in nanopores are still elusive and need to be further explored.

In this work, the adsorption and diffusion mechanisms of heavy oil four-fractions in quartz nanopores have been investigated in detail via MD simulations. Firstly, the density profiles and simulation snapshots of heavy oil fractions in quartz nanopores are used to investigate the adsorption behaviors. Furthermore, the radial distribution functions are calculated to analyze the agglomeration behaviors of heavy oil molecules. In addition, the mean square displacements evaluate the movement and diffusion of heavy oil molecules. Eventually, to further clarify the diffusion and adsorption behaviors, the self-diffusion coefficients of heavy oil molecules and the interaction energies between heavy oil fractions and quartz nanopore at different temperatures are compared to understand the desorption of heavy oil in the nanopore. This work is helpful for understanding the adsorption and desorption mechanisms of heavy oil in the nanopore.

2. Methodology

2.1. Model Construction

Heavy oil is a complex mixture with numerous hydrocarbon molecules, which can be characterized based on the broad fractions of saturates, aromatics, resins, and asphaltenes (SARA). However, the exact chemical structures of the majority of the heavy oil molecules are still unclear, and the composition and chemistry of these hydrocarbon molecules depend strongly on the geological sources and reservoir conditions. Nevertheless, some representative compounds of the four-fractions of heavy oil still can well describe the properties of most heavy oil. In this work, the eicosane (C₂₀H₄₂, MW: 282) molecule was chosen as a substitute for saturate fraction, as shown in Figure 1a. The aromatic model (C₄₆H₅₀S, MW: 635) [18] and resin model (C₄₉H₇₈S, MW: 699) [19] were selected to represent the aromatic and resin fractions, as shown in Figure 1b,c, respectively. Moreover, the asphaltenes can be structurally classified into Archipelago and Island (continental) types. The Archipelago-type model comprises several aromatic sheets attached through alkyl chains, while the Island architecture is a centered condensed aromatic molecule. In particular, the Island architecture allows the formation of aggregates by the stacking of their aromatic regions [20] and has been widely used as a basis for molecular and phenomenological models [21,22,23]. Therefore, an Island asphaltene model (C₅₄H₆₅NO₂S, MW: 792) [19] was chosen for our simulations, which contained heteroatoms of sulfur, nitrogen, and oxygen (Figure 1d). The experimentally measured contents of the SARA fractions of heavy oil in the Bohai Sea are 41.45%, 22.56%, 30.78%, and 5.21%, respectively [24]. According to the experimental data, the simulated values for the contents of the SARA fractions are 41.29%, 22.49%, 30.96%, and 5.26%, respectively. The maximum error of the mass fraction values between experiment and simulation is within ±1.0%. In addition, a classical alpha-quartz (α-SiO₂) model is applied to construct the loose sandstone nanopore in the Bohai Sea [25]. In Figure 1f, the quartz nanopore model of dimension L_x = 6.0 nm and L_y = 6.0 nm was constructed from the orthogonalized α-SiO₂ crystal (Figure 1e). According to the experimental study [26], it has been shown that the nanopores ranged from 3.0 nm to 1700.0 nm in nanometer pores of tight oil sandstone. To mimic the adsorption mechanisms at the critical nanopore size, we set the size of nanopore to be 3.0 nm in this work, as shown in Figure 1f. The nanopore surface is obtained by cleaving the (0 0 1) Miller plane of the α-SiO₂ crystal. To simulate the quartz minerals in the actual formation, the nanopore surface is hydroxylated with 9.52 hydroxyl (-OH) per nm², which falls in the range of the experimental values of SiO₂ surface (5.9~18.8 nm⁻²) [27]. To construct an reasonable adsorption model of heavy oil in a nanopore, 66 saturates, 16 aromatics, 20 resins, and 3 asphaltenes was randomly mixed in the nanopore of size 6 × 6 × 3 nm by the PACKMOL [28] and MOLTEMPLATE [29] packages, as shown in Figure 1f.

2.2. Molecular Dynamics Simulation

2.2.1. Force Fields Employed in the Simulations

In this work, the Optimized Potentials for Liquid Simulations-All Atoms (OPLS-AA) force field [30] was utilized to describe the interactions among atoms in the heavy oil phase, and the force field parameters were assigned by the LigParGen [31]. The OPLS-AA force field has been successfully utilized to simulate the adsorption properties of petroleum [11,32,33,34]. The total interaction energy of the system includes both bonded and non-bonded interactions, which can be expressed as [35]:

$$\begin{array}{l}{U}_{\mathit{bonded}}={\displaystyle \sum _{\mathit{bonds}}}{K}_r{\left(r-r_{eq}\right)}^2+{\displaystyle \sum _{\mathit{angles}}}{K}_{\theta }{\left(\theta -{\theta }_{eq}\right)}^2+\hfill \\ {\displaystyle \sum _i}\left[\frac{V_1^i}{2}\left(1+\cos \left({\phi }_i+f_{i}1\right)\right)+\frac{V_2^i}{2}\left(1-\cos \left(2{\phi }_i+f_{i}2\right)\right)+\frac{V_3^i}{2}\left(1+\cos \left(3{\phi }_i+f_{i}3\right)\right)\right]\hfill \end{array}$$

(1)

$$U_{non\text{-}bonded}=\sum _i\sum _{j>i}4{\epsilon }_{ij}\left[{\left(\frac{{\sigma }_{ij}}{r_{ij}}\right)}^{12}-{\left(\frac{{\sigma }_{ij}}{r_{ij}}\right)}^6\right]+\sum _i\sum _{j>i}\frac{q_i{q}_j}{4\pi {ϵ}_o{r}_{ij}}$$

(2)

where U_bonded and U_non-bonded are the bonded and non-bonded energies, respectively. The bonded interactions are used to model bond stretching, angle bending, and torsional bending. In Equation (1), the K_r, K_θ, and V_n are the force field constants, r and r_eq represent the length of the bond after stretching and equilibration, respectively. θ and θ_eq are the bond angle after bending and equilibration, respectively. Moreover, the third term of Equation (1) is for the torsional energy, where φ_i is the dihedral angle, V₁, V₂, and V₃ are the coefficients in the Fourier series, and f₁, f₂, and f₃ are phase angles. In addition, the non-bonded potentials consist of 12-6 Lennard-Jones (LJ) and Coulomb potentials. Herein, the LJ potential represents the van der Waals interactions, and the LJ potential parameters (ε_ij and σ_ij are the energy and size parameters between atoms i and j) between different atom types can be obtained by using the geometric mixing rule (${\epsilon }_{ij}=\sqrt{{\epsilon }_{ii}{\epsilon }_{jj}}$ and ${\sigma }_{ij}=\sqrt{{\sigma }_{ii}{\sigma }_{jj}}$). q_i is the charge of atom i, ε_o is the vacuum permittivity, and r_ij is the distance between atoms i and j. A cut-off distance of 1.4 nm was applied for the short-range forces, while the long-range Coulomb electrostatic interactions were computed using the particle-particle-particle-mesh (PPPM) method [36]. In addition, the classical CLAYFF [37] force field was utilized to describe the interaction force between hydroxylated quartz surface atoms, as shown in

Table 1. Nonbond potential parameters for CLAYFF force [37]. q represents the charge of an atom, and ε and σ represent the energy and distance between atoms.

Species	Charge, q (e)	Energy, ε (kcal/mol)	Distance, σ (Å)
H(hydroxyl)	0.4250	0.0000	0.0000
O(hydroxyl)	−0.9500	0.1554	3.5532
O(bridging)	−1.0500	0.1554	3.5532
Si	2.1000	1.8405 × 10−6	3.7950

. The CLAYFF [37] force field has been used to successfully simulate the adsorption properties of heavy oil molecules at oil—rock interactions in previous work [11,13,38].

2.2.2. MD Simulation Details

In the current work, all of the MD simulations were implemented by large-scale atomic/molecular massively parallel simulator (LAMMPS) package [39]. The periodic boundary conditions were adopted in the x and y directions of the simulation system. The time step was set to 1 femtosecond (fs). In Figure 1g, the workflow of MD simulation is as follows: firstly, the geometry optimization was carried out by the conjugate gradient method, which eliminated the interatomic overlapping. Then, an annealing process was performed to relax the interatomic residual stress, and obtain a reasonable structure under a Canonical (constant number of atoms, constant-volume, constant-temperature, NVT) ensemble with a Nosé Hoover thermostat [40]. Then, the NVT ensemble was performed for 1.0 ns to reach a specified temperature following the annealing. Subsequently, an isothermal-isobaric (constant number of atoms, constant-pressure, constant-temperature, NPT) ensemble with a Nosé Hoover thermostat and barostat [40] was conducted for 1.0 ns at the required temperature and pressure. Eventually, the MD simulation is performed 10 ns to record the required data after equilibrating the system. The time average is used instead of the ensemble average to analyze the equilibrium properties of the system.

2.3. Model Validation

To verify the stability of the simulations, the energy and density of the heavy oil system with the simulation time were plotted in Figure 2 during the NPT relaxation at different temperatures. In Figure 2a, it can be observed that the energies of the heavy oil system are stabilized around a certain value at different temperatures, and the energy gradually increases with the temperature due to the enhanced kinetic energy. In addition, the density values tend to be stable with the simulation time, but the time point at which the density value converges is different at 298 K, 398 K, and 498 K. Clearly, the density converges faster at high temperatures, while the density value is lower. As shown in Figure 2b, the converged densities of a heavy oil droplet are 0.89 g/mL, 0.84 g/mL, and 0.76 g/mL at 298 K, 398 K, and 498 K, respectively. Our density results are slightly lower than the experiment values. For instance, the experiment value of the heavy oil density in the Bohai Sea is 0.9663 g/mL at 293 K, while the simulated value is 0.89 g/mL at 298 K. In comparison with the experimental value, the error of the simulated density value is 7.9%. Despite the nuance caused by the temperature difference, the force field is an important factor affecting the accuracy of the simulated density values. Yao et al. [41] indicated that the density values obtained with different force fields are all lower than the experimental values. It is difficult to achieve the same density value as the experimental data due to the unavoidable free volume generated by the intermolecular interactions, and the differences in the physicochemical morphology between the model molecules and the real heavy oil molecules [41]. Nonetheless, these representative simulation models can still well investigate the various physical and chemical properties of heavy oil, such as adsorption [8,13,17], interface [23,42], and emulsification [43].

To further verify the simulations, the relation between pressure and density has been investigated in Figure 3. Figure 3a shows that the pressure and temperature have reached a steady state. As shown in Figure 3b, the variations of density are different with the simulation time at 2 MPa, 5 MPa, 10 MPa, and 15 MPa. It can be found that high pressure can quickly obtain a convergent density value. Moreover, a slight increase in density values with pressure can be found in Figure 3b, which is similar to the experimental results [44] (detailed in Supplementary Information Table S1). To sum up, our MD model and the applicability of the chosen calculation method are reasonable.

3. Results and Discussion

3.1. Density Profiles

As shown in Figure 4, the mass density profiles have been calculated to investigate the adsorption behaviors of heavy oil molecules in quartz nanopore at different temperatures (298 K, 398 K, and 498 K), and the density formula can be expressed as:

$${\rho }_m(z)=〈{\displaystyle \sum }_{i=1}^N\frac{M{W}_i\delta \left(z-z_i\right)}{L_x{L}_y\mathsf{\Delta }z}〉$$

(3)

where MW_i, z_i and ∆z are the weight of the heavy oil molecules, the z-coordinate of atom i and the layer thickness, respectively. The angle brackets represent the time average. In Figure 4, it can be observed that the multilayer adsorption layers have been formed near the nanopore surface. Similar phenomena had been reported by other research groups [45,46,47], such as octane in quartz nanopores [45]. As shown in Figure 4a, lots of heavy oil molecules are adsorbed on the nanopore surface at low temperature (298 K). In particular, the saturates, aromatics, and resins contribute to the density of the first adsorption layer in Figure 4a, while the asphaltene molecules are wrapped in the center of nanopore by other fractions. The weak diffusivity and strong interaction between molecules at low temperatures are the main reasons for the wrapped asphaltenes. Moreover, the peak value of the adsorption layer decreases when the temperature reaches 398 K in Figure 4b. The high temperature weakens the interactions between heavy oil molecules and increases the diffusivity of molecules, resulting in the desorption of some heavy oil molecules. It can be observed that the desorbed heavy oil molecules are mainly saturates, while the adsorption densities of the other three fractions tend to increase slightly. Furthermore, this trend can be found at 498 K in Figure 4c. In particular, it can be found that the adsorption density of aromatic molecules is most affected by temperature among the three adsorbed heavy oil fractions (aromatics, resins, and asphaltenes).

To visualize the adsorption behaviors of heavy oil molecules, the simulation snapshots of final configurations have been plotted in Figure 5 and Figure 6. From Figure 5a, it can be observed that the saturate molecules are filling the nanopore with the increase of temperature, which implies the desorption of saturate molecules. In addition, an adsorption phenomenon of aromatic molecules appears with the increase in temperature in Figure 5b. For instance, most of the aromatic molecules are wrapped by other fractions at the center of the nanopore at 298 K, while many aromatic molecules have been adsorbed on the lower surface of the nanopore at 398 K. Then, some aromatic molecules have been adsorbed on the upper and lower surfaces of the nanopore at 498 K, which indicates that the high temperature can enhance the adsorption strength of aromatic molecules. From Figure 6a, it can be observed that the resin molecules form some great resin droplets in the nanopore at 298 K. Then, the large resin droplets disappear at 398 K, and resin molecules fill the nanopore. However, some smaller resin droplets have been formed again at 498 K. The results suggest that high temperatures can reduce the interaction between resin molecules, but a higher temperature may make the agglomeration of the aromatic ring chains and alkyl chains in resin molecules (as shown in Figure 1c). In addition, the inconspicuous adsorption of asphaltene molecules has been presented with the increase in temperature in Figure 6b, which may be attributed to the large molecular weight and the wrapping of other fractions.

3.2. Radial Distribution Functions (RDFs)

To further understand the behaviors of heavy oil molecules in the nanopore, we compared the RDFs of the SARA in the heavy oil droplets and nanopore. The RDFs can be calculated by [23]:

$$g_{ab}(r)=\frac{V}{N_a{N}_b}\left({\displaystyle \sum _{i=1}^{N_a}\frac{n_{i}b(r)}{4\pi r^2\mathrm{\Delta }r}}\right)$$

(4)

where N_a and N_b represent the total numbers of atoms a and b, respectively, V stands for the volume of the simulation box, and n_ib(r) denotes the number of atom b at the radial distance of r from atom a. The characteristic carbon (C) atoms in heavy oil fractions are selected to represent the heavy oil molecules. From Equation (4), the RDFs indicate the ratio of the local density to the overall density, which represent the probability of the particles appearing at a certain position. The peak values and positions of the RDFs represent the agglomeration and distance between heavy oil molecules, respectively. As shown in Figure 7, the RDFs of the saturates in the heavy oil droplet and nanopore have been plotted at 298 K, 398 K, and 498 K. It can be observed that the peak position of the RDFs of saturates in the heavy oil droplet shift slightly to the right with the increase of temperature in Figure 7a, which indicates interactions between saturates are weakening due to high temperatures. The peak values of the RDFs of saturates exhibit a different phenomenon, that is, the peak value increases at 398 K and then decreases at 498 K. This phenomenon may be related to the curling of the chain structure at relatively high temperatures. In comparison with the saturates in the heavy oil droplet, a significant difference can be found in the temperature effect on the RDFs of saturates in the nanopore, as shown in Figure 7b. The RDFs value of saturates in nanopore significantly decreases when the temperature increases from 298 K to 398 K or 498 K. The peak position shifts slightly to the right with the increase in temperature, which is similar to the RDFs of saturates in the heavy oil droplet. In addition, we plotted the snapshots of saturate molecules in Figure 7c,d to understand the difference between the RDFs in the heavy oil droplet and nanopore at 398 K. Some large saturate clusters have been formed in the heavy oil droplet at 398 K, while the saturates still maintain a chain structure in the nanopore. The results suggest that the confined space of nanopore may suppress the agglomeration behaviors of saturates with chain structure.

Furthermore, we compared the RDFs of aromatic molecule pairs in the heavy oil droplet and the nanopore, as shown in Figure 8a,b. It can be observed that the influence of temperatures on the RDFs of aromatics in the heavy oil droplet and nanopore is the opposite. To be specific, the peak value of the RDFs of aromatics in the heavy oil droplet gradually increases with the temperature in Figure 8a, which suggests that the high temperature can induce the agglomeration of the aromatics. Figure 8c exhibits the simulation snapshots of the aromatic molecules in the heavy oil droplet at 298 K, 398 K, and 498 K, which demonstrates the agglomeration behaviors. These clusters caused by the agglomeration of aromatics gradually grow with the increase in temperature. By contrast, the RDFs values of aromatics in nanopore reduce with the increase in temperature in Figure 8b. In addition, the simulation snapshots in Figure 5b show that the aromatics in nanopore exhibit a weakened agglomeration behavior with increasing the temperature due to the adsorption effect of the nanopore surface. These results once again demonstrate that the confined space of nanopore can suppress the agglomeration behavior of aromatics induced by high temperatures.

In addition, the RDFs of the resin molecule pairs in the heavy oil droplet and the nanopore at different temperatures have been plotted in Figure 9, respectively. It is found that the effect of temperatures on the RDFs of resin molecule pairs in the heavy oil droplet and nanopore is similar. That is, the peak value of the RDFs decreases with the increase in temperature, as shown in Figure 9a,b. Furthermore, the RDFs have a significant peak shift when the temperature increases from 298 K to 398 K or 498 K, while the peak position is almost unchanged when the temperature increases from 398 K to 498 K. The results suggest that the high temperature can reduce the interactions between resin molecules, whether it is in the heavy oil droplet or nanopore. As shown in Figure 9c,d, it can be observed that the agglomeration phenomenon of resin molecules in both the heavy oil droplet and nanopore gradually disappears with the increase in temperature. However, it is found that the confined space of nanopore cannot suppress the agglomeration behaviors of resins, which may be attributed to too many chain structures (such as aromatic ring and alkane chains in Figure 1c).

Figure 10 depicts the RDFs of the asphaltene molecule pairs in the heavy oil droplet and nanopore. It can be observed that the peak value of the asphaltene molecule pairs in both the heavy oil droplet and nanopore at 398 K is higher than that of 298 K. It means that the asphaltene molecule pairs at a relatively high temperature are more likely to form the stacking structures. Similar results are also found in other research [17]. Furthermore, the peak value of the RDFs of the asphaltenes in the heavy oil droplet (Figure 10a) decreases and the peak position shifts to right at 498 K, and this phenomenon is more obvious in the RDFs of asphaltene pairs in nanopore (Figure 10b). According to the corresponding simulation snapshots, some stacking types of asphaltene molecules in the heavy oil droplet and nanopore have been found at different temperatures. In Figure 10a, a side-to-face stacking and a large distance between asphaltenes (in the red rectangle) were found in the heavy oil droplet at 298 K, which result in the formation of three low peaks. Furthermore, the face-to-face stacking (in the green circle) was formed at 398 K. However, the side-to-side stacking (in the blue circle) between asphaltene molecules appears when the temperature increases from 398 K to 498 K, which suggests that the asphaltene molecule pairs at a quite high temperature are likely to break the face-to-face stacking. In addition, as shown in Figure 10b, the asphaltene molecule pairs in the nanopore form a clear side-to-face stacking (in the red rectangle), which may be attributed to the constraint of the nanopore. Furthermore, asphaltene molecules in the nanopore form a twined side-to-face stacking (in the green circle) due to the confined space at 398 K. In comparison with the face-to-face stacking of asphaltenes in the heavy oil droplet, it implies that the nanopore can suppress the agglomeration behaviors of asphaltenes induced by high temperature. Furthermore, the scattered asphaltene molecules (in the blue circle) in the nanopore again verify the effect of the confined space of the nanopore on the agglomeration behavior of asphaltene molecules.

3.3. Mean Square Displacement (MSD)

The mobility of heavy oil molecules can be evaluated by the MSD functions calculated as follows [21]:

$$\mathrm{MSD}\left(t\right)=〈{\left|{\mathit{r}}_i(t)-{\mathit{r}}_i(0)\right|}^2〉$$

(5)

where r_i(t) is the position of molecule i at time t and r_i(0) is the initial position of molecule i. The angle brackets indicate the time average. The MSD is described as the square of the difference between the position of a molecule at time t and the initial position, providing a good understanding of the behaviors of the molecules in terms of displacements [23]. As shown in Figure 11 and Figure 12, the MSD curves of saturate, aromatic, resin, and asphaltene molecules in the heavy oil droplet and nanopore were presented at different temperatures, respectively. It can be observed that the movements of heavy oil molecules in oil droplet are an order of magnitude higher than that of heavy oil molecules in the nanopore. In addition, the movements of saturate, aromatic, resin, and asphaltene molecules in the heavy oil droplet increase with temperature in Figure 11. The effect of temperatures on the movements of heavy oil molecules in nanopore has a slight decrease compared with the oil molecules in the heavy oil droplet. These results imply that the mobility of heavy oil molecules in a nanopore may be suppressed due to the confined space and adsorption effect.

3.4. Self-Diffusion Coefficient and Adsorption Energy

The diffusivity of heavy oil molecules in nanopore is closely related to the adsorption effect, which influences the detachment of heavy oil from nanopore. Herein, we have analyzed the relations between the self-diffusion coefficients (D_i) and interaction energy. The D_i can be calculated by fitting the MSD based on the Einstein formula [21,48]:

$$D_i=\frac{1}{6}\underset{t\to \mathit{\infty }}{\lim }\frac{d}{dt}\sum _{i=1}^N〈{\left|{\mathit{r}}_i^c(t)-{\mathit{r}}_i^c(0)\right|}^2〉$$

(6)

where N represents the total number of particles. Anomalous diffusion of heavy oil molecules might occur with non-negligible influence. Therefore, according to the Einstein formula, the trajectories of MD simulations were extracted over the time frame when the Einstein diffusion occurs (that is, the slope of log (MSD) vs. log (Time) is 1) to calculate the D_i [21]. In addition, the interaction energy (E_oil/quartz) of the SARA fractions is calculated to evaluate the interaction strength by the following equation [49]:

$$E_{oil/quartz}=\frac{E_{total}-(E_{oil}+E_{quartz})}{N_{oil}}$$

(7)

where E_oil/quartz indicates the interaction energy between oil molecules and quartz. E_oil and E_quartz represent the energy of heavy oil fraction and quartz nanopore, respectively. The E_total is the total energy of the whole system containing heavy oil fraction and quartz nanopore, and the N_oil represents the number of heavy oil molecules. The negative interaction energy value means there is an attraction between heavy oil and quartz nanopore.

As shown in Figure 13a, the D_i values of the SARA fractions in nanopore have been calculated at 298K, 398 K, and 498 K. It can be found that the D_i values of saturates, aromatics, and resins increase with temperature, while the D_i of the asphaltenes decreases and then increases from 298 K to 498 K. The results indicate that high temperatures can enhance the diffusivity of most heavy oil molecules in the nanopore. Nevertheless, the diffusivity of asphaltenes at 398 K is slightly suppressed, which may be mainly attributed to the twined side-to-face stacking of asphaltenes in nanopore at 398 K, as shown in Figure 10b. Furthermore, to explore the relations between diffusivity and adsorption, the corresponding E_oil/quartz of the SARA fractions have been calculated by Equation (7) in Figure 13b. It can be observed that the E_oil/quartz of saturates decrease with the increase in temperature, this trend is opposite to its D_i. The result suggests that high temperatures can reduce the adsorption of saturates and increase their diffusivity, and the variations of D_i and E_oil/quartz weaken from 398 K to 498 K. This phenomenon can also be found in the simulation snapshots in Figure 5a and the RDFs in Figure 7b. Unlike saturates, the E_oil/quartz and D_i of aromatics and resins both increase with temperature in Figure 13. The competitive phenomenon between adsorption and diffusivity of aromatics and resins results in their different adsorption behaviors. For instance, the aromatic molecules in nanopore exhibit a competitive adsorption behavior between upper and lower rock from 298 K to 498 K in Figure 5b. In addition, the resin molecules in nanopore exhibit a similar agglomeration behavior to that in the heavy oil droplet due to the competitive behavior between adsorption and diffusivity, as shown in Figure 9c,d. Moreover, the E_oil/quartz of asphaltenes show a non-monotonic trend with temperature, which is similar to the D_i. That is, the adsorption strength of asphaltenes enhances and then weakens from 298 K to 398 K. The results suggest that the adsorption of asphaltenes is stronger than diffusivity at 398 K due to the twined side-to-face stacking, while the adsorption and diffusivity of asphaltenes remain a competitive behavior at 298 K and 498 K.

In addition, the E_oil/quartz in Figure 13b indicates the average E_oil/quartz of SARA molecules, which exhibits a comprehensive effect including the adsorption between SARA fractions and the adsorption of nanopore surface. To only consider the intrinsic adsorption of a single heavy oil fraction, we have calculated the E_oil/quartz of single saturate, aromatic, resin, and asphaltene molecules in nanopore at 298 K, 398 K, and 498 K respectively. As shown in Figure 14, it can be observed that the adsorption strength of single asphaltene molecule is relatively stronger than other fractions. In comparison with Figure 13b, the E_oil/quartz of the asphaltene molecule decreases continuously with increasing temperature due to the disappeared stacking. On the contrary, the E_oil/quartz of the resin molecule gradually increases with temperature, which considers the stacking between the flexible aromatic ring chains of a single resin molecule (Figure 1c). Nevertheless, the non-monotonic E_oil/quartz of the aromatic molecule may be limited by its inflexible aromatic ring chain (Figure 1b). In addition, the E_oil/quartz of a single saturate molecule is insensitive to temperatures, which may be attributed to its nonpolar property and light structure. It is suggested that the adsorption strength of single heavy oil fractions is closely related to their molecular structures.

4. Conclusions

In this paper, we have investigated the influence of temperature on the adsorption and diffusion behaviors of heavy oil fractions via MD simulations. The major conclusions are summarized as follows:

(1)

Heavy oil forms a denser multilayer adsorption layer on the nanopore surface, and temperature can alter the adsorption behaviors of heavy oil fractions. As the temperature increases, the saturate molecules will be desorbed from the nanopore surface, the aromatic molecules will be aggregated in the adsorption layer, the resin molecules will remain a similar aggregation behavior in the nanopore, and the asphaltene molecules will aggregate towards the nanopore surface.

(2)

The confined space of the nanopore will suppress the agglomeration behaviors of the heavy oil molecules. The saturates in nanopore maintain a chain structure, while the saturates in the heavy oil droplet forms some clusters. The aromatics also maintain a similar phenomenon to saturates. Nevertheless, the confined space of nanopore cannot suppress the agglomeration behaviors of resin molecules due to too many chain structures. The agglomeration of resin molecules in both the heavy oil droplet and nanopore gradually disappears with increasing temperatures. The face-to-face stacking of asphaltene molecules in the heavy oil droplet is changed to the side-to-face (twined) stacking in nanopore suggesting that the confined space can decrease the stacking of asphaltene molecules.

(3)

The influence of temperature on the movements of heavy oil molecules in a nanopore has a slight decrease compared with the oil molecules in the heavy oil droplet, which suggests that the mobility of heavy oil molecules in a nanopore may be suppressed due to the confined space and adsorption effect.

(4)

The adsorption and diffusivity of aromatic, resin, and asphaltene molecules in nanopore exhibit a competitive phenomenon with the increase in temperature, resulting in different adsorption behaviors of heavy oil molecules.

This work provides molecular-level theoretical insights into the understanding of the adsorption and diffusion mechanisms of heavy oil fractions in quartz nanopores. Further work will focus on the simulation of the desorption and displacement mechanisms of heavy oil fractions with strong adsorption in nanopores during EOR. This will be beneficial to improve the efficient exploitation of remaining oil adsorbed in nanopores.