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.