Surface Properties Govern Oil Transport in Nanochannels: Molecular Dynamics Insights

Background

With global energy demand rising and conventional resources depleting, unconventional shale oil has attracted intense research interest because of its vast reserves [1,2]. Despite the fact that shale‑oil resources exceed conventional oil by more than threefold [2], the exploitable fraction remains limited due to its confinement within nanoscopic pores ranging from 2 to 100 nm [4,5]. This confinement induces pronounced surface effects that alter fluid properties—such as density, wettability, and diffusion—leading to transport behaviours that differ markedly from macroscopic systems [10–12]. For example, water flows faster in carbon nanotubes than in equivalent macroscale tubes [6,7], and the hydrophilic/hydrophobic character of carbon nanotubes shifts with diameter [9]. Prior molecular dynamics work has shown that shear stress and pore size influence water and oil dynamics in carbon and silica nanoslits [13–15]. However, systematic studies on how surface chemistry and topography govern shale‑oil transport are scarce. Understanding these effects is crucial for optimizing extraction technologies and has broader implications for membrane separations, biomolecular transport, and energy storage devices [16–20].

Methods

All simulations were performed using the Discover module of Materials Studio (Accelrys Inc.) with the COMPASS all‑atom force field, which accurately captures bonded, cross‑term, and non‑bonded interactions in condensed phases:

$$ {E}_{\mathrm{total}}={E}_{\mathrm{valence}}+{E}_{\mathrm{cross}\text{-}term}+{E}_{\mathrm{nonbond}} \tag{1} $$

Full expressions for each term are given in the original publication. Van der Waals forces were modeled with a Lennard‑Jones potential (cut‑off 15.5 Å), and electrostatics were treated via the Ewald method. Simulations were run in the NVT ensemble at 298 K using the Andersen thermostat. Periodic boundaries were applied in all directions. Data were sampled every 5 ps.

The silica substrate, representing the dominant mineral in shale, was constructed from the (001) surface and hydroxylated to mimic geologic conditions [26]. Two opposing surfaces, each 1.5 × 7 × 0.85 nm³, created a slit channel with widths of 2, 4, and 6 nm. Octadecane, pyridine, and phenol—common shale‑oil constituents—were introduced to achieve a bulk density of 0.8 g cm⁻³ (40 molecules in the 2‑nm channel; 407, 344, 80, and 120 molecules for pyridine, phenol, and octadecane in the 4‑ and 6‑nm channels, respectively). After energy minimization and a 500‑ps equilibration, a constant body force equivalent to 3.1 × 10⁻¹⁴ N per atom was applied along the channel axis to drive the molecules, a standard approach in MD transport studies [27–29].

Results and Discussion

Effect of Channel Width

Under the imposed force, the number of molecules crossing the channel increased steadily (Supporting Information, Fig. S1). Figure 2 demonstrates that wider channels yield larger COM displacements after 2 ns. The COM position is defined as

$$ {z}_{\mathrm{COM}}={\displaystyle \sum_i\frac{m_i}{M}{r}_i} \tag{5} $$

In the 2‑nm channel, the COM moved only 0.85 nm, whereas in the 6‑nm channel the displacement exceeded 30 nm. This pronounced difference arises from stronger adsorption in narrower pores, as quantified by the average interaction energy

$$ {E}_{\mathrm{average}\ \mathrm{interaction}}=\frac{E_{\mathrm{total}}-(E_{\mathrm{oil}}+E_{\mathrm{substrate}})}{N} \tag{6} $$

Figure 3 shows a clear inverse correlation: higher interaction energy leads to reduced displacement, with a threefold increase in energy causing a >30‑fold drop in COM shift. Interaction effects diminish as the channel widens, underscoring the size‑dependent nature of nanoscale transport.

Layered density profiles (Fig. 2) reveal a 5‑Å thick interfacial layer adjacent to the walls. Octadecane orientation within these layers is predominantly parallel to the surface (θ ≈ 80°–90°), while molecules in the channel centre exhibit random orientations. Diffusion coefficients, derived from mean‑square displacements (Eq. 7), vary with radial position: they are highest at the centre and decrease towards the wall (Fig. 5a). The 6‑nm channel displays a steep parabolic profile, whereas the 2‑nm channel shows nearly uniform, low diffusion, resembling a piston‑like front. These trends can be attributed to the competition between wall adsorption (which suppresses motion in the near‑surface “shadow” region) and inter‑molecular momentum transfer.

Effect of Polarity

Simulations in a 4‑nm silica channel with comparable numbers of phenol, pyridine, and octadecane illustrate the role of molecular polarity. Despite identical forces, octadecane’s COM displacement is ~16 × greater than that of phenol or pyridine (Fig. 8). Total interaction energies (vdW + electrostatic) are markedly higher for the polar species, with electrostatics contributing roughly half of the total (Fig. 8). Dipole moments calculated via first‑principle methods confirm that octadecane (0.032 Debye) is virtually non‑polar, whereas phenol (1.306 Debye) and pyridine (2.245 Debye) possess significant dipoles. Consequently, polar molecules experience stronger binding to the silica surface, reducing mobility.

Effect of Material Composition

Comparisons among silica, gold, and calcite channels (Fig. 9) show that gold strongly immobilises octadecane, whereas calcite allows moderate displacement. Figure 10 indicates that the COM shift correlates with the average interaction energy: silica exhibits the lowest energy and highest displacement, gold the opposite, and calcite intermediate. The similar energies between silica and calcite yet differing displacements suggest that surface atom type and texture also influence transport.

Effect of Surface Roughness

Introducing nanoscale cavities (3 Å or 6 Å deep) into the silica surface alters local flow. Cavities host trapped oil molecules, which reduce the velocity of adjacent fluid layers (Fig. 11). COM displacements for 3‑Å and 6‑Å cavities are 3.95 Å and 3.07 Å, respectively—slightly higher than the flat‑surface value (3.17 Å) for the shallow cavity, but lower for the deeper cavity. The enhancement with shallow cavities likely results from an effective widening of the channel, while deeper cavities impose stronger adsorption and hinder flow.

Conclusions

This molecular‑dynamics study demonstrates that nano‑channel width, oil polarity, substrate material, and surface roughness collectively dictate oil transport. The 6‑nm channel affords a >30‑fold greater COM displacement than the 2‑nm channel, and the centre‑of‑channel diffusion coefficient nearly doubles that near the wall. Polar molecules and high‑interaction substrates markedly suppress mobility, while shallow surface cavities can modestly enhance transport. These insights are essential for designing nano‑structured extraction systems and for broader applications where fluid transport in confined geometries is critical.

Abbreviations

COM:

Center of mass

COMPASS:

Condensed‑phase optimized molecular potential for atomistic simulation studies

MD:

Molecular dynamics

vdW:

van der Waals