Fig. 3 Snapshots showing the initial configurations of the
seven P-H slit pore models of 5 nm width with varying water
concentration. Cw stands for the water concentration. The clay structure
uses the same color codes as in Fig. 1 . The color codes of
fluid molecules are: ethane in light green; dodecane in light blue; and
H2O in red. In the first panel, the direction of the
arrows indicates the direction of an imposed acceleration in subsequent
sections of this paper.
2.2 Discussion of the Force
Field
We apply the ClayFF force field to
describe the interatomic interactions for illite and the
cations53 as is common with clay molecular
simulations54,55. Water molecules are described using
a flexible SPC model and the shake algorithm is used to make the two O-H
bonds and the H-O-H angle rigid56. We use the OPLS
All-Atom force field to represent the organic components (ethane and
dodecane)57. We also use the Lorentz-Berthelot mixing
rules to calculate the interaction parameters between unlike atoms. It
should be noted that the mixing rules of ClayFF and OPLS All-Atom force
fields are different. The OPLS All-Atom force field uses a geometric
mixing rule to obtain the Lennard-Jones interaction parameters between
unlike atoms57, while the Lorentz-Berthelot mixing
rule is used for ClayFF force field53. For the
interaction parameters between clay and organic components, the use of
Lorentz-Berthelot mixing rule has been documented in Hao et
al.50, Wu et al.58, and Zhao et
al.59 for methane adsorption in clay minerals and is
also applied in our study.
2.3 Simulation Details
We use the Large-scale
Atomic/Molecular Massively Parallel Simulator (LAMMPS)60 with periodic boundary conditions in all three
directions with short-range interactions represented by a Lennard-Jones
(LJ) 12-6 term.61 The cutoff distance for the
short-range nonbonded van der Waals interactions is 8 Å and the
long-range electrostatic interactions are calculated by the
Fourier-based Ewald summation method62 - the
particle-particle/particle-mesh (PPPM) method with a precision value of
10−6 63.
Our workflow is as follows: In the initial set-up, hydrocarbon and water
are placed randomly within the pores using the Packmol
package64. We then run equilibrium MD (EMD)
simulations using an NPT ensemble using a time step of 0.01 fs for 100ps
and of 1 fs for the next 5 ns. Pressure is controlled at 400 atm by the
Parrinello-Rahman barostat65 while temperature is held
at 350 K by the Nose Hoover thermostat 66. With
equilibrium in the NPT ensemble, we switch to the NVT ensemble and
continue the simulations for another 5 ns.
After EMD simulations, non-equilibrium MD (NEMD) simulations are
performed for another 10 ns to mimic hydrocarbon-water transport.
Several methods exist for inducing flow in molecular dynamics, such as
forced67–69, surface-induced, pressure
difference70, osmotic pressure71 and
gravitational approaches72–74. In our simulations, we
adopt the gravitational technique because other methods have been known
to cause a distortion of the velocity-field and density profile along
the flow direction74–76. A constant external
gravity-like force, parallel to the basal plane (along the x-direction
as shown in Fig. 3), is exerted on all hydrocarbon-water molecules
inside the illite nanopores77,78. To ensure a linear
system response, we maintain a constant acceleration of each particle on
the order of 10-4 to 10-3nm/ps2 79,80. In our study, the range of our
acceleration is from 0.0005 to 0.002 nm/ps2.
It should be noted that the fluid temperature in MD simulations is
normally calculated from the kinetic energy via summation of the
velocity-squared of all particles in the system81.
However, the particle velocities along the direction of driving force
are made up of two components: thermal velocity and center-of-mass
velocity (the imposed streaming velocity). To ensure that the
center-of-mass velocity does not contribute to the fluid temperature,
only the velocity components perpendicular to the driving force are used
to calculate the fluid temperature 76. The temperature
and pressure as a function of simulation time during NEMD simulations
are presented in Fig. S2 in the Supporting Information.