Fig. 3 Initial configurations of the seven PH slit pore models of 5 nm width with varying water concentration. Cw stands for the water concentration. The color codes of clay are as in Fig. 1. In the first panel, the direction of the arrows indicates the direction of an imposed acceleration in subsequent sections of this paper. The fluid-color regimes: ethane (light green), H2O (red) and dodecane (light blue).

2.2 Simulation Details

Our modeling relies on the use of the ClayFF and OPLS All-Atom force field to describe illite53,54 and organic components55respectively, while the flexible SPC model and the shake algorithm are applied to model water. Lorentz-Berthelot mixing rules describe the interactions between different atoms. We use LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) 56 with periodic boundary conditions applied in 3 directions. A more detailed description is available in our previous work42,48.
Our workflow is as follows: In the initial set-up, fluid molecules are placed randomly between the slit pore using the Packmol package57. We then run equilibrium MD (EMD) simulations using an NPT ensemble for 5 ns. After equilibrium of isothermal–isobaric (NPT) ensemble, we further run the simulations for another 5 ns with a canonical (NVT) ensemble. Temperature and pressure are maintained at 400 atm and 350 K using Parrinello-Rahman barostat58 and Nose Hoover thermostat59 respectively.
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 forced60–62, surface-induced, pressure difference63, osmotic pressure64 and gravitational approaches65–67. 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 direction67–69. A uniform external force (along the x-direction as shown in Fig. 3), is imposed on all molecules inside the illite nanopores70,71. Published literature suggests using small accelerations on the order of 10-4 to 10-3 nm/ps2 72,73 to avoid velocity jumps in the system and therefore, in our study, the acceleration ranges from 0.0005 to 0.002 nm/ps2.
It should be noted that in MD simulations, fluid temperature is calculated from the kinetic energy74. However, in NEMD simulations, the molecule velocities in the direction of the driving force consists of both thermal velocities and the imposed center-of-mass velocity. To ensure the temperature is maintained, our fluid temperature calculations do not include the center-of-mass velocity69. The temperature and pressure during the NEMD simulations are shown in Fig. S2 in the Supporting Information.