*2.1. Modeling Approach*

Molecular dynamics rely on decoupling between the nucleus and electron motions. While the nucleus motion is classically treated through Newton's mechanics, the e ffects of the electron motions are considered through partial charges placed on the molecular structure. The force field is a set of equations that describe the molecular interactions and includes intermolecular and intramolecular interactions. Intermolecular interactions describe the interactions among atoms from di fferent molecular entities and include both the Van der Waals (VDW) and electrostatics interactions. While the Lenard Jones potential (Equation (1)) is one of the widely-used models to model the VDW interactions, electrostatic interactions are modeled by Coulomb's law (Equation (2)). On the other hand, the intramolecular interactions maintain the molecular entity and include the bond, angle, dihedral and improper interactions. Given the initial configuration of the system and the molecular interactions,

$$\mu(r)\_{Lf} = 4\varepsilon \left[ \left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r} \right)^{6} \right] \tag{1}$$

$$\left(u(r)\_{\text{coulomb}} = \frac{Q\_1 Q\_2}{4\pi \varepsilon\_0 r}\right) \tag{2}$$

where ε is the depth of the potential well (dimensions: ML2T−2), σ is the distance at which the potential energy is zero (dimension: L), *Q* is the atomic charge (dimension: AT) and 0 is the permittivity of free space (dimensions: <sup>M</sup>−1L−3T4A2).

## *2.2. Simulation Details*

We used graphite to represent the pore material, decane and pentane to represent the hydrocarbons and carbon dioxide, respectively, and nitrogen and methane to represent the EOR's gases. We elected to simulate rough solid surfaces, which realistically capture the complexity of the porous media. As shown in Figure 1, our system consists of a silt pore connected to a reservoir containing the EOR's gas. This reservoir is constrained by a wall that acts as a fixed boundary or constant-rate boundary. In addition, a production region is defined on the other side of the pore. The whole simulation matrix is shown in Table 1. The dimensions of the simulation cell are lx = 504 Å, ly = 42.6 Å and lz = 90/70 Å. The solid substrate comprise three sheets of graphite. We used equal molar gas density for all gases. Note that equal gas densities do not correspond to equal fluid pressures, as shown in Figure 2.

**Figure 1.** System setup and the simulation evolution: (**a**) the initial setup of the system where the hydrocarbons are placed inside the pore and the Enhanced Oil Recovery (EOR) gas (in this case CO2) is placed outside; (**b**) the hydrocarbon equilibration; (**c**) the EOR gas equilibration; (**d**) the system extended to include a production region, where the pore was opened for production and a constant-velocity wall is allowed; (**e**) the system after 0.1 ns; (**f**) the system after 0.5 ns; (**g**) the system after 1 ns. Color code: red is graphite, blue is gas and yellow is hydrocarbon.


**Table 1.** Force field parameters of the system components.

**Figure 2.** The impact of pressure on the bulk density of carbon dioxide, nitrogen and methane at temperature = 350 K (National Institute of Standards and Technology (NIST) data [51]).

We used the OPLS\_UA force field [52] to model the hydrocarbon interactions and united-atom force fields for the EOR's gases as detailed in Table 2. We used Moltemplate software to set up the initial system [53] and LAMMPS to run all the simulations [54]. We started by annealing the hydrocarbons inside the pore for 0.1 ns using canonical ensemble and then equilibrating the EOR's gas for 0.1 ns. After that, the simulation was extended to include a production region (an extra 250 Å in the x-direction). In this scenario, we simulated concurrent displacement, where both the injected fluid and hydrocarbons are moving in the same directions. For constant-rate boundary cases, we applied a constant velocity of 0.0002 Kcal/mole-Angstrom on the wall. We ran the simulation for 1 ns and recorded the hydrocarbon production. We created rough surfaces by carving out grooves inside the solid substrate. We used the Lorentz–Berthelot mixing rules to model the cross interactions. We used the Nosé–Hoover thermostat with a damping parameter of 100 time steps.


**Table 2.** Simulation matrix of the scenarios studied.

We also simulated the Huff and Puff process, where the EOR's gas is kept in contact with the hydrocarbon system for 10 ns. After that, the region of the EOR's gas is evacuated to allow the hydrocarbon recovery. This scenario depicts countercurrent displacement. We recorded the production for another 5 ns.
