**2. Results and Discussion**

### *2.1. Emulsified Crude Oil Droplet*

It was believed that SDS increased the hydrophilicity of oil droplets by increasing the hydrophilic surface area of the droplets [34]. In order to study the surface condition of oil droplets with a different SDS content in each system, the solvent-accessible surface area (SASA) was calculated and shown in Figure 1. With an increasing number of SDS molecules, we noted that both the hydrophilic surface and the hydrophobic surface of crude oil droplets also increased. Meanwhile, the ratio of hydrophilic area to hydrophobic area also increased significantly. Therefore, adding SDS molecules could increase the hydrophilic surface area of oil droplets more significantly. In addition, the greater the number of SDS molecules, the greater the hydrophilic surface area. oil droplets also increased. Meanwhile, the ratio of hydrophilic area to hydrophobic area also increased significantly. Therefore, adding SDS molecules could increase the hydrophilic surface area of oil droplets more significantly. In addition, the greater the number of SDS molecules, the greater the hydrophilic surface area.

droplets with a different SDS content in each system, the solvent-accessible surface area (SASA) was calculated and shown in Figure 1. With an increasing number of SDS molecules, we noted that both the hydrophilic surface and the hydrophobic surface of crude

*Molecules* **2022**, *27*, x FOR PEER REVIEW 3 of 15

**Figure 1.** Average hydrophobic and hydrophilic surface areas of emulsified oil droplet. **Figure 1.** Average hydrophobic and hydrophilic surface areas of emulsified oil droplet.

### *2.2. Dynamic Behavior of Oil Droplets under BPEF 2.2. Dynamic Behavior of Oil Droplets under BPEF*

In order to study the behavior of oil droplets under electric field, BPEF with *E* = 0.50 V/nm was applied in the z-direction of all systems. Figure 2 displayed the conformational changes of oil droplets with differing SDS content during the electric field output stage. As can be seen from Figure 2, all oil droplets gradually deformed under the electric field, elongated in the z-direction and migrated toward the opposite direction of the electric field. Moreover, SDS and asphaltene molecules concentrated at the end of the oil droplet. Excess SDS molecules were distributed along the entire surface of the deformed oil droplets (System IV, V). To clearly see the distribution of SDS and asphaltene, oil droplets of System II were partly zoomed in. It can be seen that the SDS and asphaltene molecules aggregated at the head of the oil droplet's moving direction, with the negative sulfonic acid groups of SDS and carboxyl groups of asphaltene molecules facing the opposite direction of the electric field. Therefore, we thought it was the polar SDS and asphaltene In order to study the behavior of oil droplets under electric field, BPEF with *E* = 0.50 V/nm was applied in the z-direction of all systems. Figure 2 displayed the conformational changes of oil droplets with differing SDS content during the electric field output stage. As can be seen from Figure 2, all oil droplets gradually deformed under the electric field, elongated in the z-direction and migrated toward the opposite direction of the electric field. Moreover, SDS and asphaltene molecules concentrated at the end of the oil droplet. Excess SDS molecules were distributed along the entire surface of the deformed oil droplets (System IV, V). To clearly see the distribution of SDS and asphaltene, oil droplets of System II were partly zoomed in. It can be seen that the SDS and asphaltene molecules aggregated at the head of the oil droplet's moving direction, with the negative sulfonic acid groups of SDS and carboxyl groups of asphaltene molecules facing the opposite direction of the electric field. Therefore, we thought it was the polar SDS and asphaltene molecules that guided the movement of oil droplets under the electric field.

molecules that guided the movement of oil droplets under the electric field. In Figure 2 we also noted that the states of two oil droplets in five systems were different at 400 ps. In System II and V, the two oil droplets collided at 400 ps. Whereas, for System I, III and IV, this didn't occur. To investigate the impact of SDS concentration on the coalescence of oil droplets driven by electric field, the collision time was summarized in Figure 3. A collision occurred when the minimum distance between two oil droplets was less than 0.35 nm. It was found that the addition of SDS molecules can reduce the In Figure 2 we also noted that the states of two oil droplets in five systems were different at 400 ps. In System II and V, the two oil droplets collided at 400 ps. Whereas, for System I, III and IV, this didn't occur. To investigate the impact of SDS concentration on the coalescence of oil droplets driven by electric field, the collision time was summarized in Figure 3. A collision occurred when the minimum distance between two oil droplets was less than 0.35 nm. It was found that the addition of SDS molecules can reduce the collision time of oil droplets, especially with the 6.2% SDS concentration condition

collision time of oil droplets, especially with the 6.2% SDS concentration condition.

*Molecules* **2022**, *27*, x FOR PEER REVIEW 4 of 15

**Figure 2.** The behavior of oil droplets under BPEF output duration. Asphaltene molecules were colored pink. SDS molecules were colored yellow. Resins and light crude oil molecules were colored green. **Figure 2.** The behavior of oil droplets under BPEF output duration. Asphaltene molecules were colored pink. SDS molecules were colored yellow. Resins and light crude oil molecules were colored green. **Figure 2.** The behavior of oil droplets under BPEF output duration. Asphaltene molecules were colored pink. SDS molecules were colored yellow. Resins and light crude oil molecules were colored green.

**Figure 3.** Collision time of oil droplets in five systems under the electric field of 0.5 V/nm. **Figure 3.** Collision time of oil droplets in five systems under the electric field of 0.5 V/nm.

### *2.3. Surface Charge Distribution*

*2.3. Surface Charge Distribution*  The electrostatic potential surface of the oil droplet can reflect its charge redistribution under electric field. Considering the two oil droplets in each system are the same, The electrostatic potential surface of the oil droplet can reflect its charge redistribution under electric field. Considering the two oil droplets in each system are the same, only

only one droplet's electrostatic potential was calculated. The electrostatic potential dia-

**Figure 3.** Collision time of oil droplets in five systems under the electric field of 0.5 V/nm.

The electrostatic potential surface of the oil droplet can reflect its charge redistribution under electric field. Considering the two oil droplets in each system are the same, only one droplet's electrostatic potential was calculated. The electrostatic potential dia-

one droplet's electrostatic potential was calculated. The electrostatic potential diagrams were obtained for different systems at the initial and specific time during the simulation (Figure 4). It can be found that under the influence of the hydrophilic and negatively charged asphaltene molecules on the surface, some areas of the oil droplets appear electronegative (blue area) before electric field is applied. The electronegative area increases with the increase of SDS content. However, the electrostatic potential at the surface of the deformed oil droplet noticeably changed under the electric field, which is manifested as one end of ellipsoidal oil droplet being electronegative toward the opposite direction of the electric field and the other end being electropositive, as in System I and II. These revealed that the redistribution of the charge of oil droplets under electric field resulted in the droplets' polarization. This phenomenon was consistent with the experimental observation that the charge of the oil droplets under the electric field is positive in the direction of the electric field, and negative in the opposite direction. Meanwhile, we found that for Systems III, IV and V, the oil droplets' polarization was not obvious. To explain this, the number density of SDS in oil droplets under electric field was analyzed at the same time. In Figure 5, we defined the middle of the oil drop as 0 and the moving direction as the positive direction. It can be seen that with an increase in SDS, it tended to distribute on the surface of the whole deformed oil droplet, which could further explain why the electronegativity area of the deformed oil droplet increased with the increase in SDS.

The dynamic behavior of SDS and asphaltene molecules of oil droplets and the electrostatic potential distribution on the surface of oil droplets displayed that the mobile negative charges on oil droplets moved toward the opposite direction of the applied electric field. However, what causes the two oil droplets moving in the same direction to collide, and its relationship with SDS content is unclear. Figure 6 presented the centroid distance between two oil droplets and the average elongation length *le* of the two oil droplets along the *z* direction from the application of the electric field to the collision of oil droplets in each system. We can find that even when the two oil droplets were deformed under electric field, the centroid distance between the two oil droplets remained approximately 10 nm in all systems. This means that due to the oil droplets having the same composition in each system, they moved along the opposite direction of the electric field at almost the same speed, so they kept almost the same initial centroid distance. However, the average elongation length *le* of the two oil droplets in the five systems studied were significantly different. It was found that in all systems, the oil droplets start length was about 6 nm in diameter, and their length increased with time; the average elongation length *le* of the oil droplets exceeded 10 nm near the collision time point. This means that when the length of the oil droplet is stretched enough, the two oil droplets are connected head to end; that is, a collision occurs. Meanwhile, we noted that in Figure 6b, the order of the growth rate of the average elongation length *le2* from largest to smallest is System II, System V, System III ≈ System IV and System I, which is similar to the trend showing the variation of the collision time of studied systems in this work. Therefore, for the O/W emulsion system with uniform distribution of oil droplets, we thought the demulsification collision time in the electric field is significantly affected by SDS. Adding appropriate SDS surfactant into O/W systems can effectively reduce the power consumption.

**Figure 4.** Electrostatic potential surface of the oil droplets during simulation. **Figure 4.** Electrostatic potential surface of the oil droplets during simulation. **Figure 4.** Electrostatic potential surface of the oil droplets during simulation.

**Figure 5.** Density profile of SDS in oil droplets under electric field at the same time. **Figure 5.** Density profile of SDS in oil droplets under electric field at the same time. **Figure 5.** Density profile of SDS in oil droplets under electric field at the same time.

**Figure 6.** Centroid distance between two oil droplets (**a**) and average elongation length *le* of the two oil droplets along the z direction (**b**). **Figure 6.** Centroid distance between two oil droplets (**a**) and average elongation length *le* of the two oil droplets along the *z* direction (**b**).

As discussed above, SDS and asphaltene molecules guide the entire oil droplet to

move in the opposite direction of the electric field. It was predicted that the average elon-

asphaltene molecules. Thus, we calculated the MSD of SDS and asphaltene molecules for the five systems in Figure 7. It was found that the order of SDS and asphaltene molecules' diffusion from largest to smallest was System II, System V, System IV, System III and System IV in the five systems studied; this was consistent with the order of the average elongation length of oil droplets under electric field. We thought that in System I the asphaltene molecules interacted more strongly with the surrounding oil molecules due to the

As discussed above, SDS and asphaltene molecules guide the entire oil droplet to move in the opposite direction of the electric field. It was predicted that the average elongation length of oil droplets in electric field is related to the diffusivity of the SDS and asphaltene molecules. Thus, we calculated the MSD of SDS and asphaltene molecules for the five systems in Figure 7. It was found that the order of SDS and asphaltene molecules' diffusion from largest to smallest was System II, System V, System IV, System III and System IV in the five systems studied; this was consistent with the order of the average elongation length of oil droplets under electric field. We thought that in System I the asphaltene molecules interacted more strongly with the surrounding oil molecules due to the influence of its structure, which decreased its mobility under electric field. However, the negative SDS molecules are smaller and demonstrate strong mobility in the electric field, thus increasing their overall mobility. However, this does not mean that the greater the SDS content in the oil droplets, the greater the mobility of negatively charged molecules. Therefore, the SDS content of the oil droplets have great significance on the demulsification effect. Meanwhile, we calculated the root-mean-square fluctuation (RMSF) of oil droplets during the electric field output stage (as shown in Supplementary Materials: Figure S1). By comparing the RMSF of the three systems, we found that the fluctuation of System II and System IV was stronger than that of System I. The addition of SDS could have accelerated the movement of oil droplets, which was similar to the calculation result of MSD. *Molecules* **2022**, *27*, x FOR PEER REVIEW 8 of 15 influence of its structure, which decreased its mobility under electric field. However, the negative SDS molecules are smaller and demonstrate strong mobility in the electric field, thus increasing their overall mobility. However, this does not mean that the greater the SDS content in the oil droplets, the greater the mobility of negatively charged molecules. Therefore, the SDS content of the oil droplets have great significance on the demulsification effect. Meanwhile, we calculated the root-mean-square fluctuation (RMSF) of oil droplets during the electric field output stage (as shown in Supplementary Materials: Figure S1). By comparing the RMSF of the three systems, we found that the fluctuation of System II and System IV was stronger than that of System I. The addition of SDS could have accelerated the movement of oil droplets, which was similar to the calculation result of MSD.

**Figure 7.** Mean square displacements of SDS and asphaltene molecules in five systems. **Figure 7.** Mean square displacements of SDS and asphaltene molecules in five systems.

### *2.4. Aggregation Behavior of Oil Droplets*

*2.4. Aggregation Behavior of Oil Droplets*  The purpose of electric field demulsification is to aggregate dispersed oil droplets to achieve oil/water separation. Conformations of the oil droplet at the beginning of the collision were selected as the initial structure to simulate the behavior of oil droplets after the shut-off of the electric field (see Figure 8). It can be seen that the oil droplets in contact with each other can continue to aggregate even in the absence of electric field. Taking System II as an example, we found that some asphaltene and SDS molecules, which guided the movement of the oil droplets, formed a contact surface between the two oil The purpose of electric field demulsification is to aggregate dispersed oil droplets to achieve oil/water separation. Conformations of the oil droplet at the beginning of the collision were selected as the initial structure to simulate the behavior of oil droplets after the shut-off of the electric field (see Figure 8). It can be seen that the oil droplets in contact with each other can continue to aggregate even in the absence of electric field. Taking System II as an example, we found that some asphaltene and SDS molecules, which guided the movement of the oil droplets, formed a contact surface between the two oil droplets

droplets and then migrated to the surface of the oil droplets under the influence of hydro-

during the aggregation of oil droplets in the five systems. (as shown in Figure S2). We found that the radius of gyration of the oil droplets gradually decreased. Therefore, the

droplets that collided would gradually aggregate into a whole.

and then migrated to the surface of the oil droplets under the influence of hydrophilic groups. At the same time, the hydrophobic components inside the interfacing oil droplets aggregated into a whole. Meanwhile, we calculated the radius of gyration (Rg) during the aggregation of oil droplets in the five systems. (as shown in Figure S2). We found that the radius of gyration of the oil droplets gradually decreased. Therefore, the droplets that collided would gradually aggregate into a whole. *Molecules* **2022**, *27*, x FOR PEER REVIEW 9 of 15

**Figure 8.** The change in the behaviors of oil droplets during the electric field shut-off period in five systems. **Figure 8.** The change in the behaviors of oil droplets during the electric field shut-off period in five systems.

#### *2.5. Mechanism of Aggregation of Oil Droplets 2.5. Mechanism of Aggregation of Oil Droplets*

It had been proposed [35] that the surface charges of oil droplets will rearrange under BPEF. The positive and negative charges at the adjacent areas of the two oil droplets are opposite under the action of the electric field. Thus, the adjacent areas of oil droplets always attract each other along the BPEF direction. Here, we verified and explained the accumulation mechanism of oil droplets through theoretical methods. We calculated the interaction energy of the two oil droplets in all systems with *E* = 0.50 V/nm during the whole process from dispersion to aggregation. The calculated results are shown in Figure 9. The potential energy of the interaction between the two oil droplets is divided into two parts. The cyan areas represented the change in the potential energy between the oil droplets from dispersion to collision (i.e., electric field output durations), and the blue areas represented the change in the potential energy with time from collision to aggregation (i.e., electric field shut-off durations). At the same time, we calculated the root-meansquare deviation (RMSD) of crude oil droplets in the five systems (Figure S3). It can be seen from the figure that the aggregated oil droplets were basically stable after 4.0 ns. We found that the potential energy of the electrostatic interaction between the oil droplets was almost 0 kJ/mol during the entire electric field application process. The potential energy of the van der Waals interactions between the two oil droplets from dispersion to collision was also almost 0 kJ/mol in the output electric field stage. However, the potential energy of the van der Waals interactions between the oil droplets noticeably decreased after collision during the aggregation process (electric field shut-off stage). It means that in the electric field demulsification process, the adjacent areas of the two oil droplets with opposite charges have no obvious effect on the attraction and aggregation of the oil droplets. The van der Waals forces between the oil droplets are the main force in the demulsi-It had been proposed [35] that the surface charges of oil droplets will rearrange under BPEF. The positive and negative charges at the adjacent areas of the two oil droplets are opposite under the action of the electric field. Thus, the adjacent areas of oil droplets always attract each other along the BPEF direction. Here, we verified and explained the accumulation mechanism of oil droplets through theoretical methods. We calculated the interaction energy of the two oil droplets in all systems with *E* = 0.50 V/nm during the whole process from dispersion to aggregation. The calculated results are shown in Figure 9. The potential energy of the interaction between the two oil droplets is divided into two parts. The cyan areas represented the change in the potential energy between the oil droplets from dispersion to collision (i.e., electric field output durations), and the blue areas represented the change in the potential energy with time from collision to aggregation (i.e., electric field shut-off durations). At the same time, we calculated the root-mean-square deviation (RMSD) of crude oil droplets in the five systems (Figure S3). It can be seen from the figure that the aggregated oil droplets were basically stable after 4.0 ns. We found that the potential energy of the electrostatic interaction between the oil droplets was almost 0 kJ/mol during the entire electric field application process. The potential energy of the van der Waals interactions between the two oil droplets from dispersion to collision was also almost 0 kJ/mol in the output electric field stage. However, the potential energy of the van der Waals interactions between the oil droplets noticeably decreased after collision during the aggregation process (electric field shut-off stage). It means that in the electric field demulsification process, the adjacent areas of the two oil droplets with opposite charges have no obvious effect on the attraction and aggregation of the oil droplets. The van der Waals forces between the oil droplets are the main force in the demulsification process.

fication process.

**Figure 9.** Interaction energy and its decomposed composition between two oil droplets during simulation. **Figure 9.** Interaction energy and its decomposed composition between two oil droplets during simulation.

### **3. Methods and Materials**

### **3. Methods and Materials**  *3.1. Simulation Details*

*3.1. Simulation Details*  All MD simulations were performed in GROMACS 2019.6 software package. The GROMOS 53a6 force field [36] was used. The force field parameters of oil droplet composition were generated by the Automated Topology Builder (ATB) [37,38]. The simple point charge (SPC) model was selected for water molecules. The parameters of sodium ions All MD simulations were performed in GROMACS 2019.6 software package. The GROMOS 53a6 force field [36] was used. The force field parameters of oil droplet composition were generated by the Automated Topology Builder (ATB) [37,38]. The simple point charge (SPC) model was selected for water molecules. The parameters of sodium ions (Na<sup>+</sup> ) that neutralize negative charge have been discussed in the literature [39].

(Na+) that neutralize negative charge have been discussed in the literature [39]. Each system was energy-minimized using the steepest descent method before the simulation. The NVT ensemble at 300 K was performed with velocity rescaling thermostat. The NPT ensemble at 0.1 MPa and 300 K was performed with Berendsen pressure coupling. In the simulation, velocity rescaling thermostat with a time constant of 0.1 ps was selected as the temperature coupling method, and Berendsen pressure coupling with a time constant of 1.0 ps was selected as the pressure coupling method; the isothermal compression factor was set to 4.5 × 10−5 bar−1. The periodic boundary condition was applied along three dimensions. During the simulation, van der Waals interactions used Lennard−Jones 12-6 potential, and the cutoff was set to 1.4 nm. The Coulombic interaction used particle-mesh Ewald (PME) summation method. The initial velocities were assigned according to Maxwell−Boltzmann distribution. The time step chose 2 fs. The trajectory Each system was energy-minimized using the steepest descent method before the simulation. The NVT ensemble at 300 K was performed with velocity rescaling thermostat. The NPT ensemble at 0.1 MPa and 300 K was performed with Berendsen pressure coupling. In the simulation, velocity rescaling thermostat with a time constant of 0.1 ps was selected as the temperature coupling method, and Berendsen pressure coupling with a time constant of 1.0 ps was selected as the pressure coupling method; the isothermal compression factor was set to 4.5 <sup>×</sup> <sup>10</sup>−<sup>5</sup> bar−<sup>1</sup> . The periodic boundary condition was applied along three dimensions. During the simulation, van der Waals interactions used Lennard−Jones 12-6 potential, and the cutoff was set to 1.4 nm. The Coulombic interaction used particlemesh Ewald (PME) summation method. The initial velocities were assigned according to Maxwell−Boltzmann distribution. The time step chose 2 fs. The trajectory was saved every 10 ps. VMD 1.9.3 was used for trajectory visualization.

### was saved every 10 ps. VMD 1.9.3 was used for trajectory visualization. *3.2. Simulation Systems*

### *3.2. Simulation Systems*  3.2.1. Molecular Models of Crude Oil

3.2.1. Molecular Models of Crude Oil Owing to the high complexity of crude oil, especially for asphaltene and resins, the asphaltenes demonstrate a key role in the stabilization of water-in-crude oil emulsions and significantly impact the rheological properties of crude oil [40]. Two types of asphaltene (i.e., the number of each type of asphaltene is four) and six types of resins [41] (i.e., the number of each type of resin is five) were selected based on previous studies, as shown in Figure 10. In addition to asphaltenes and resin molecules, four types of alkanes (32 Owing to the high complexity of crude oil, especially for asphaltene and resins, the asphaltenes demonstrate a key role in the stabilization of water-in-crude oil emulsions and significantly impact the rheological properties of crude oil [40]. Two types of asphaltene (i.e., the number of each type of asphaltene is four) and six types of resins [41] (i.e., the number of each type of resin is five) were selected based on previous studies, as shown in Figure 10. In addition to asphaltenes and resin molecules, four types of alkanes (32 hexane, 29 heptane, 34 octane and 40 nonane molecules), two types of cyclanes (22 cyclohexane and 35 cycloheptane molecules) and two types of aromatics (13 benzene and 35 toluene molecules) were

selected as light oil components, referring to Song and Miranda's work [41–43]. Moreover, the concentration of resins and asphaltene in the crude oil was about 38%, which met the content of heavy oil components in crude oil [41]. toluene molecules) were selected as light oil components, referring to Song and Miranda's work [41–43]. Moreover, the concentration of resins and asphaltene in the crude oil was about 38%, which met the content of heavy oil components in crude oil [41].

hexane, 29 heptane, 34 octane and 40 nonane molecules), two types of cyclanes (22 cyclohexane and 35 cycloheptane molecules) and two types of aromatics (13 benzene and 35

*Molecules* **2022**, *27*, x FOR PEER REVIEW 11 of 15

**Figure 10.** Asphaltene and resin molecules used in simulation. **Figure 10.** Asphaltene and resin molecules used in simulation.

### 3.2.2. Emulsified Oil Droplet

3.2.2. Emulsified Oil Droplet First, the components of crude oil including alkanes, cyclanes, aromatics, asphaltenes and resins were randomly inserted into a cubic box (x = 10 nm, y = 10 nm, z = 10 nm). To eliminate overlapping, energy minimization was then performed. After that, a 30 ns NPT First, the components of crude oil including alkanes, cyclanes, aromatics, asphaltenes and resins were randomly inserted into a cubic box (x = 10 nm, y = 10 nm, z = 10 nm). To eliminate overlapping, energy minimization was then performed. After that, a 30 ns NPT ensemble simulation was performed to obtain a reasonable density. The equilibrium configuration after NPT run was shown in Figure 11a.

ensemble simulation was performed to obtain a reasonable density. The equilibrium configuration after NPT run was shown in Figure 11a. Second, the above crude oil was then solvated in an 8 nm × 8 nm × 8 nm simulation Second, the above crude oil was then solvated in an 8 nm × 8 nm × 8 nm simulation box with 19,230 water molecules. Energy minimization and a 20 ns NVT ensemble MD simulation was carried out to obtain the emulsified oil droplet (Figure 11b).

box with 19,230 water molecules. Energy minimization and a 20 ns NVT ensemble MD simulation was carried out to obtain the emulsified oil droplet (Figure 11b). Third, emulsified oil droplets with different amounts of SDS adsorbed on their surface was constructed. SDS micelles were constructed using Packmol. The above spherical oil droplets were then placed in the center of a new box (10 nm × 10 nm × 15 nm) and SDS Third, emulsified oil droplets with different amounts of SDS adsorbed on their surface was constructed. SDS micelles were constructed using Packmol. The above spherical oil droplets were then placed in the center of a new box (10 nm × 10 nm × 15 nm) and SDS micelles were placed close to oil droplets (Figure 11c). Then, Na<sup>+</sup> counter ions and solvent were added. After energy minimization and a 20 ns NVT simulation, emulsified oil droplet systems were derived (Figure 11d).

micelles were placed close to oil droplets (Figure 11c). Then, Na+ counter ions and solvent were added. After energy minimization and a 20 ns NVT simulation, emulsified oil drop-

let systems were derived (Figure 11d).

**Figure 11.** Schematic diagram of the building of emulsified oil droplet. (**a**) Equilibrium configuration of bulk crude oil after NPT simulation, (**b**) emulsified oil droplet in water, (**c**) initial configuration of emulsified oil droplet and SDS micelle, (**d**) Equilibrium configuration of emulsified oil droplet with SDS. For clarity, water molecules in (**c**,**d**) were not shown. (**e**) Lateral view of model simulation for the movement and aggregation behavior of oil droplets in O/W emulsion under BPEF. For clarity, water molecules in the system were not shown. **Figure 11.** Schematic diagram of the building of emulsified oil droplet. (**a**) Equilibrium configuration of bulk crude oil after NPT simulation, (**b**) emulsified oil droplet in water, (**c**) initial configuration of emulsified oil droplet and SDS micelle, (**d**) Equilibrium configuration of emulsified oil droplet with SDS. For clarity, water molecules in (**c**,**d**) were not shown. (**e**) Lateral view of model simulation for the movement and aggregation behavior of oil droplets in O/W emulsion under BPEF. For clarity, water molecules in the system were not shown.

We assumed that the oil droplets distributed in the emulsion had the following conditions: First, the centroids of the two oil droplets were approximately along the *z*-axis direction; Second, the centroids of the two drops were about 10.0 nm apart; Finally, two identical emulsified oil drop models with counter ions were placed in a 10 × 10 × 50 nm3 box with a separation distance of about 10 nm, as shown in Figure 11e. Afterward, water molecules were added to solvate the system. The energy minimization and a 10 ns NVT simulation were applied to ensure emulsion system equilibrium. Subsequently, BPEF was imposed on all systems to study coalescence of the two droplets. The composition of each emulsified oil droplet system was shown in Table 1. We assumed that the oil droplets distributed in the emulsion had the following conditions: First, the centroids of the two oil droplets were approximately along the *z*-axis direction; Second, the centroids of the two drops were about 10.0 nm apart; Finally, two identical emulsified oil drop models with counter ions were placed in a 10 <sup>×</sup> <sup>10</sup> <sup>×</sup> 50 nm<sup>3</sup> box with a separation distance of about 10 nm, as shown in Figure 11e. Afterward, water molecules were added to solvate the system. The energy minimization and a 10 ns NVT simulation were applied to ensure emulsion system equilibrium. Subsequently, BPEF was imposed on all systems to study coalescence of the two droplets. The composition of each emulsified oil droplet system was shown in Table 1.


**Table 1.** Details of the emulsified oil droplet systems. **Table 1.** Details of the emulsified oil droplet systems.

### V 1 60 68 46,680 28.5% **4. Conclusions**

**4. Conclusions**  In this paper, molecular dynamics simulations were performed to study the behavior of oil droplets in O/W emulsion. The differences in oil droplets emulsified by different amounts of SDS were compared. Three major conclusions were derived. First, the hydro-In this paper, molecular dynamics simulations were performed to study the behavior of oil droplets in O/W emulsion. The differences in oil droplets emulsified by different amounts of SDS were compared. Three major conclusions were derived. First, the hydrophilicity of oil droplets increases with increasing SDS content in the oil droplet. When

philicity of oil droplets increases with increasing SDS content in the oil droplet. When

electric field is applied, oil droplets move in the opposite direction of the electric field. The molecules in the oil droplet underwent redistribution. SDS and asphaltene with negatively charged functional groups were transferred to the head of the droplet along the direction of movement. The electrostatic potential surface of the oil droplet proved that the BPEF made the molecules redistribute in the droplet, which resulted in its surface potential redistribution as well. This is consistent with the theoretical hypothesis proposed by this experiment. Meanwhile, the collision time of oil droplets in all simulation systems was different due to the different SDS mass fraction, and the collision time was the shortest for the oil droplets with 6.2% SDS. The average elongation length *le* of the two oil droplets along the *z* direction explained that SDS molecules could change the elongation length of the oil droplets in the electric field. The MSD of SDS and asphaltene molecules under electric field showed that the mobility was the strongest in System II. Therefore, the elongation length of the oil droplets in System II was the largest, and this system was the least time consuming. Second, the oil droplets after collision can self-aggregate after electric field shut-off. SDS and asphaltene molecules on the contact surface between the two oil droplets migrated to the surface of the oil droplets under the influence of hydrophilic groups. Lastly, the adjacent areas of the two oil droplets with opposite charges have no obvious effect on the attraction and aggregation of the oil droplets, and the van der Waals forces between oil droplets are the main force in the demulsification process.

**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/molecules27082559/s1, Figure S1: Root-mean-square fluctuation (RMSF) of oil droplets in system I, system II and system IV, Figure S2: Radius of gyration (Rg) of oil droplets in five systems, Figure S3: Root-mean-square deviation (RMSD) of oil droplets in five systems.

**Author Contributions:** S.L. did the molecular dynamic simulations and wrote the manuscript. S.Y. supervised the molecular dynamic simulations and revised the manuscript. H.Z. designed and carried out the MD simulations and revised the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Natural Science Foundation of Shandong Province (No. ZR2021MB040).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** We gratefully appreciate the financial support from the Natural Science Foundation of Shandong Province.

**Conflicts of Interest:** The authors declare no conflict of interest.

**Sample Availability:** Not applicable.

### **References**


**Ding Li 1,2, Shuixiang Xie <sup>1</sup> , Xiangliang Li <sup>3</sup> , Yinghua Zhang <sup>3</sup> , Heng Zhang 2,\* and Shiling Yuan <sup>2</sup>**


**Abstract:** CO<sup>2</sup> enhanced oil recovery (CO<sup>2</sup> -EOR) has become significantly crucial to the petroleum industry, in particular, CO<sup>2</sup> miscible flooding can greatly improve the efficiency of EOR. Minimum miscibility pressure (MMP) is a vital factor affecting CO<sup>2</sup> flooding, which determines the yield and economic benefit of oil recovery. Therefore, it is important to predict this property for a successful field development plan. In this study, a novel model based on molecular dynamics to determine MMP was developed. The model characterized a miscible state by calculating the ratio of CO<sup>2</sup> and crude oil atoms that pass through the initial interface. The whole process was not affected by other external objective factors. We compared our model with several famous empirical correlations, and obtained satisfactory results—the relative errors were 8.53% and 13.71% for the two equations derived from our model. Furthermore, we found the MMPs predicted by different reference materials (i.e., CO2/crude oil) were approximately linear (R<sup>2</sup> = 0.955). We also confirmed the linear relationship between MMP and reservoir temperature (TR). The correlation coefficient was about 0.15 MPa/K in the present study.

**Keywords:** minimum miscible pressure; CO<sup>2</sup> enhanced oil recovery; molecular dynamics

### **1. Introduction**

Global warming has caused great changes such as continued sea level rise, which is irreversible over hundreds to thousands of years. CO<sup>2</sup> is the culprit of this phenomenon. CCUS (CO<sup>2</sup> capture, utilization, and storage) is a new technology developed from CCS (CO<sup>2</sup> capture and storage) that can bring economic benefits while reducing CO<sup>2</sup> emissions and alleviating global warming [1]. CO<sup>2</sup> enhanced oil recovery (CO2-EOR) is one of the effective ways of CCUS. The captured CO<sup>2</sup> is squeezed into the oil reservoirs that have been exploited, and the interaction between CO<sup>2</sup> and crude oil is used to improve the properties of the crude oil, thereby displacing more crude oil from the crust [2]. Research has shown that CO2-EOR can improve crude oil recovery significantly and extend the life of oil reservoirs [3,4]. Hence, CO2-EOR has been fundamentally well researched in laboratories and applied in industries as an efficient approach since the 1970s [5].

There are two different miscible and immiscible states in CO2-EOR. Under the former condition, CO<sup>2</sup> and crude oil can completely integrate into one phase, resulting in a much higher recovery rate than the latter. For the former, there is a minimum pressure above which CO<sup>2</sup> and crude oil can be miscible. This minimum pressure value, also called the minimum miscible pressure (MMP), is a vital parameter in the process of CO2-EOR. Nevertheless, considering the massive influencing factors, the accurate determination of MMP remains a major challenge [6].

**Citation:** Li, D.; Xie, S.; Li, X.; Zhang, Y.; Zhang, H.; Yuan, S. Determination of Minimum Miscibility Pressure of CO2–Oil System: A Molecular Dynamics Study. *Molecules* **2021**, *26*, 4983. https://doi.org/10.3390/ molecules26164983

Academic Editor: James Gauld

Received: 7 May 2021 Accepted: 16 August 2021 Published: 17 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

To date, there are various ways to predict MMP such as experimental measurement and computational methods. The former has been widely used due to their high precision. Within them, slim-tube experiments [7–9], as a necessary test in the industry, is considered to be the standard experimental procedure. Rising-bubble apparatus (RBA) [10,11] and vanishing interfacial tension (VIT) [12–15] are also frequently utilized to determine MMP because of their simplicity and flexibility. Although these experimental measurements have accurate techniques, they still suffer from some disadvantages including time-consumption and operation cost. Furthermore, it is difficult for any experimental method to simulate the real conditions of the crude oil reservoirs completely so that their results are greatly influenced by the instruments.

The application of computational techniques is an available alternative approach to experiments. In 1960, the first empirical MMP correlation was proposed by Benham et al. [16]. The reported equation was correlated using three pseudo-components presenting a multi-components system, and some satisfactory results were obtained based on this model. Thereafter, an increasingly number of correlations were developed for MMP prediction [17–20]. Researchers found that the more useful parameters an equation used, the better performance the model had [21]. These parameters generally included reservoir temperature (TR), composition of drive gas (CO2, H2S, N2, and C1–C5), molecular weight of C5+ fraction in crude oil (MWc5+), and the ratio of volatile (C<sup>1</sup> and N2) to intermediate (C2–C4, H2S, and CO2) in crude oil (Vol./Int.).

In addition to the conventional empirical formula models, the parameters above are often used in some intelligent algorithms based on machine learning. For instance, artificial neural networks (ANNs) can learn from large amounts of input data, and reflect their relationships more effective than conventional techniques [22]. Determination of network structure and its parameters are two crucial steps in achieving high performance from ANN. One part of the data is used to train and look for a suitable structure and optimal parameters, while the other tests the prediction accuracy of the model. Based on the principle, back propagation (BP) [23] and radial basis function (RBF) [24] are proposed. Beyond that, a series of optimization methods such as genetic algorithm (GA) [25], particle swarm optimization (PSO) [26], support vector machine (SVM) [27], and hybrid-ANFIS [28] have also been developed for MMP determination. In a previous study [29], we compared four estimation methods and found that the machine learning intelligent algorithm had a higher precision to the MMP than pure linear model. In addition, some reports that combined multiple approaches showed better results [30–34].

However, all of the above methods cannot give a direct explanation of the MMP from a microscopic view. They are all based on the existing oilfield data, which means that the established model will inevitably be affected by specific situation. To put it another way, these methods can be considered as pure mathematical statistics methods that have low levels of universality for different CO2-EOR.

Against this backdrop, the current study proposes a novel MMP prediction model at the molecular level, and the research process was not affected by other external objective factors. Therefore, the model represents a new strategy. First, we built a simulation box that contained CO<sup>2</sup> and crude oil with an obvious phase interface. To mimic the contact between CO<sup>2</sup> and crude oil, these molecules were gradually mixed until they were miscible with time evolution by using molecular dynamics. After calculating the ratio of CO2/crude oil atoms that passed through the initial interface, we found the connection between the ratio value with the miscible state. When the ratio changes from decreasing to stable, it indicates that the system has entered a miscible state, while the pressure corresponding to the inflection point is MMP. Figure 1 is a flow chart that shows all the main steps of modeling. The main objective of this study was to reveal the principle of the MMP formation at the molecular level and provide more feasible ideas for the prediction of the MMP.

**Figure 1.** Flowchart of proposed MMP prediction model. (i) Construction of the simulation system, (ii) Extracting number density data after MD simulation, (iii) Determination of initial miscible time, (iv) Reconfirmation of initial miscible time, (v) Processing data from initial miscible time to the end, (vi) Processing data from P1 to P6 at T, (vii) Acquisition of MMP. **Figure 1.** Flowchart of proposed MMP prediction model. (i) Construction of the simulation system, (ii) Extracting number density data after MD simulation, (iii) Determination of initial miscible time, (iv) Reconfirmation of initial miscible time, (v) Processing data from initial miscible time to the end, (vi) Processing data from P1 to P6 at T, (vii) Acquisition of MMP.

#### **2. Simulation Method 2. Simulation Method**

#### *2.1. Simulation and Force Field 2.1. Simulation and Force Field*

The molecular dynamics simulation was performed by the GROMACS 4.6.7 package [35,36], and AMBER 03 all-atom force field [37]. Parameters set for all components of crude oil and CO2 were generated from Automated Topology Builder and Repository databases [38,39]. The molecular dynamics simulation was performed by the GROMACS 4.6.7 package [35,36], and AMBER 03 all-atom force field [37]. Parameters set for all components of crude oil and CO<sup>2</sup> were generated from Automated Topology Builder and Repository databases [38,39].

The convergence criterion of energy minimization was 1000 kJ/(mol·nm). In the simulation, a velocity rescaling thermostat with a 0.1 ps time constant was selected as the temperature coupling method [40]. Berendsen pressure coupling with 1.0 ps time constant was selected as the pressure coupling method. The isothermal compression factor was set to 4.5 × 10−5 bar−1 [41]. The time step was 2 fs, and periodic boundary conditions were applied in the XY directions [42]. Walls were set at the top and bottom of the Z-direction in the simulated box to ensure that all atoms passed through the initial interface to achieve the miscibility. Bond lengths were constrained by the LINCS algorithm [43]. During the simulation, van der Waals interactions with the Lennard–Jones potential was cut off at 1.4 nm. Coulomb interaction used the particle-mesh Ewald summation method [44,45]. The Verlet list was updated every 10 steps. The Maxwell–Boltzmann distribution was employed to set the initial atomic velocities of the systems [46]. The trajectories were integrated by the leapfrog Verlet algorithm [47]. The convergence criterion of energy minimization was 1000 kJ/(mol·nm). In the simulation, a velocity rescaling thermostat with a 0.1 ps time constant was selected as the temperature coupling method [40]. Berendsen pressure coupling with 1.0 ps time constant was selected as the pressure coupling method. The isothermal compression factor was set to 4.5 <sup>×</sup> <sup>10</sup>−<sup>5</sup> bar−<sup>1</sup> [41]. The time step was 2 fs, and periodic boundary conditions were applied in the XY directions [42]. Walls were set at the top and bottom of the Z-direction in the simulated box to ensure that all atoms passed through the initial interface to achieve the miscibility. Bond lengths were constrained by the LINCS algorithm [43]. During the simulation, van der Waals interactions with the Lennard–Jones potential was cut off at 1.4 nm. Coulomb interaction used the particle-mesh Ewald summation method [44,45]. The Verlet list was updated every 10 steps. The Maxwell–Boltzmann distribution was employed to set the initial atomic velocities of the systems [46]. The trajectories were integrated by the leapfrog Verlet algorithm [47].

### *2.2. Simulation System*

*2.2. Simulation System*  In a real situation, the chemical components of crude oil are highly complex. Under the current experimental conditions, it is time-consuming to precisely analyze the exact constitution of its components. In order to get as close as possible to the real situation, the oil model was designed based on Miranda's works [48,49], which were used to explore the interface properties between crude oil and different fluids. Their model contained alkanes (72 hexane, 66 heptane, 78 octane, and 90 nonane molecules), cyclanes (48 cyclohexane and 78 cycloheptane molecules), and aromatics (30 benzene and 78 toluene molecules), and has been proven reliable by Song et al. [50]. In a real situation, the chemical components of crude oil are highly complex. Under the current experimental conditions, it is time-consuming to precisely analyze the exact constitution of its components. In order to get as close as possible to the real situation, the oil model was designed based on Miranda's works [48,49], which were used to explore the interface properties between crude oil and different fluids. Their model contained alkanes (72 hexane, 66 heptane, 78 octane, and 90 nonane molecules), cyclanes (48 cyclohexane and 78 cycloheptane molecules), and aromatics (30 benzene and 78 toluene molecules), and has been proven reliable by Song et al. [50].

At 333 K and 10 MPa, all alkanes, cyclanes, and aromatics were added into a cubic box (x = 9 nm, y = 9 nm, z = 9 nm) randomly. Then, energy minimization was performed to eliminate opposed-conformation. In order to mimic the state of crude oil in the reservoir, we performed a 30 ns NPT ensemble simulation to obtain its equilibrium state. After equilibration, the size of simulation box changed to 5.2 nm × 5.2 nm × 5.2 nm. box (x = 9 nm, y = 9 nm, z = 9 nm) randomly. Then, energy minimization was performed to eliminate opposed-conformation. In order to mimic the state of crude oil in the reservoir, we performed a 30 ns NPT ensemble simulation to obtain its equilibrium state. After equilibration, the size of simulation box changed to 5.2 nm × 5.2 nm × 5.2 nm. Furthermore, we built a box of the same size, stochastically adding 561 CO2 molecules to mimic the supercritical CO2 fluid (333 K, 10 MPa). Energy minimization and 30 ns NVT

At 333 K and 10 MPa, all alkanes, cyclanes, and aromatics were added into a cubic

Furthermore, we built a box of the same size, stochastically adding 561 CO<sup>2</sup> molecules to mimic the supercritical CO<sup>2</sup> fluid (333 K, 10 MPa). Energy minimization and 30 ns NVT ensemble simulation enabled the CO<sup>2</sup> to reach its equilibrium state. To simulate the contact between CO<sup>2</sup> and crude oil, the two boxes were integrated into one rectangular simulation box, and the height of new box in the Z-direction was slightly increased to 11.2 nm to avoid intermolecular overlap, as shown in Figure 2. After that, at least 10 ns NPT ensemble MD simulation was performed. ensemble simulation enabled the CO2 to reach its equilibrium state. To simulate the contact between CO2 and crude oil, the two boxes were integrated into one rectangular simulation box, and the height of new box in the Z-direction was slightly increased to 11.2 nm to avoid intermolecular overlap, as shown in Figure 2. After that, at least 10 ns NPT ensemble MD simulation was performed.

*Molecules* **2021**, *26*, x FOR PEER REVIEW 4 of 14

**Figure 2.** Construction of the simulation system. (**a**) CO<sup>2</sup> phase. (**b**) Crude oil phase. (**c**) Initial simulation box. (**d**) Integral of CO2/crude oil molecules passing through the initial phase interface.

### **Figure 2.** Construction of the simulation system. (**a**) CO2 phase. (**b**) Crude oil phase. (**c**) Initial simulation box. (**d**) Integral of CO2/crude oil molecules passing through the initial phase interface. **3. Results and Discussion**

**3. Results and Discussion**  In the last NPT ensemble simulation, the size of the box gradually stabilized over time. When CO2 was miscible with crude oil, the NPT ensemble achieved equilibrium and the box size remained unchanged. However, the change in the size of the box cannot In the last NPT ensemble simulation, the size of the box gradually stabilized over time. When CO<sup>2</sup> was miscible with crude oil, the NPT ensemble achieved equilibrium and the box size remained unchanged. However, the change in the size of the box cannot intuitively reflect the miscibility process. Therefore, we introduced the number density of CO2/crude oil.

intuitively reflect the miscibility process. Therefore, we introduced the number density of CO2/crude oil. The product of number density and volume of a box is the total number of atoms. When the system achieves equilibrium, the ratio of CO2/crude oil atoms in the upper or lower half of the box should be 50%. In our simulation system, the Z-direction height will not change due to the presence of walls. Therefore, the integral change in the number density in the Z-direction (i.e., the integral bars in Figure 2d) can reflect the change in the The product of number density and volume of a box is the total number of atoms. When the system achieves equilibrium, the ratio of CO2/crude oil atoms in the upper or lower half of the box should be 50%. In our simulation system, the Z-direction height will not change due to the presence of walls. Therefore, the integral change in the number density in the Z-direction (i.e., the integral bars in Figure 2d) can reflect the change in the size of the box and further reflect the mixing progress. When the system achieves equilibrium, the integral of the number density in Z-direction will also be constant.

### size of the box and further reflect the mixing progress. When the system achieves *3.1. Definition of Initial Miscible Time*

equilibrium, the integral of the number density in Z-direction will also be constant. *3.1. Definition of Initial Miscible Time*  First, the initial miscible time was defined. It refers to the moment when CO2 and the crude oil phases just reach the miscible state during their mixing progress, and they can keep the miscible state afterward. The purpose is to ensure that data after this time are miscibility data. We used the CO2 phase as an example to illustrate the calculations. Its number density data were extracted along the Z-direction of the box from 0.5 ns to 10.0 ns First, the initial miscible time was defined. It refers to the moment when CO<sup>2</sup> and the crude oil phases just reach the miscible state during their mixing progress, and they can keep the miscible state afterward. The purpose is to ensure that data after this time are miscibility data. We used the CO<sup>2</sup> phase as an example to illustrate the calculations. Its number density data were extracted along the Z-direction of the box from 0.5 ns to 10.0 ns every 0.5 ns after the NPT ensemble was run. Hence, there were 20 sets of data in total. We can obtain the number of CO<sup>2</sup> atoms in the lower half of the box by integrating the density of CO<sup>2</sup> along the Z-direction below the initial interface at each cut-off time. Since CO<sup>2</sup> was not distributed in the box below the initial interface at the beginning, the integral

values we obtained corresponded to the number of CO<sup>2</sup> atoms passing through the initial interface at the cut-off time. More vividly, it is an integral bar in the three-dimensional space of the box along the Z-direction, as shown in Figure 2d. molecule always kept in continuous random motion, thus it is normal to have positive and negative fluctuations after miscibility. For the selection of the initial miscibility time, the establishment standard is to find the time when the curve becomes stable and the

every 0.5 ns after the NPT ensemble was run. Hence, there were 20 sets of data in total. We can obtain the number of CO2 atoms in the lower half of the box by integrating the density of CO2 along the Z-direction below the initial interface at each cut-off time. Since CO2 was not distributed in the box below the initial interface at the beginning, the integral values we obtained corresponded to the number of CO2 atoms passing through the initial interface at the cut-off time. More vividly, it is an integral bar in the three-dimensional

Furthermore, a curve of the number of CO2 atoms passing through the initial

in the lower half of the box at 333 K and 10 MPa: it gradually increased from zero to a stable value (about 49.09), and then tended to be stable. It is worth noting that each

Furthermore, a curve of the number of CO<sup>2</sup> atoms passing through the initial interface over time can be plotted. Figure 3 shows the change in the number of CO<sup>2</sup> atoms in the lower half of the box at 333 K and 10 MPa: it gradually increased from zero to a stable value (about 49.09), and then tended to be stable. It is worth noting that each molecule always kept in continuous random motion, thus it is normal to have positive and negative fluctuations after miscibility. For the selection of the initial miscibility time, the establishment standard is to find the time when the curve becomes stable and the change is very gentle after the miscibility reference line (i.e., the 49.09 line in Figure 3a). This is also the time when the miscibility has just been achieved. The first-order variance of data with time evolution (Figure 3b) reflects the trend of data changes more intuitively. It should be noted that it is not the "initial miscibility" that is already zero, but the time corresponding to the point relatively close to zero. Based on the situation in Figures 2 and 3, it can be guaranteed that the time at 4 ns: (i) the vertical axis value is already very close to the reference line, and (ii) the curve's upward trend has slowed down. change is very gentle after the miscibility reference line (i.e., the 49.09 line in Figure 3a). This is also the time when the miscibility has just been achieved. The first-order variance of data with time evolution (Figure 3b) reflects the trend of data changes more intuitively. It should be noted that it is not the "initial miscibility" that is already zero, but the time corresponding to the point relatively close to zero. Based on the situation in Figures 2 and 3, it can be guaranteed that the time at 4 ns: (i) the vertical axis value is already very close to the reference line, and (ii) the curve's upward trend has slowed down. Based on similar treatments, crude oil atoms passing through the initial interface with time evolution (Figure S1).) and its first-order variance (Figure S1b) can also be drawn. Figures 3 and S1 show that the changes in CO2 and crude oil were quite similar. Combined with the analyses above, we can preliminarily conclude that 4 ns is the initial miscible time, and is also the time when the CO2–oil system achieved miscibility.

**Figure 3.** The number of CO2 atoms passing through the initial interface (**a**) and its first–order variance (**b**) with time evolution (333 K, 10 MPa). **Figure 3.** The number of CO<sup>2</sup> atoms passing through the initial interface (**a**) and its first–order variance (**b**) with time evolution (333 K, 10 MPa).

*3.2. Reconfirmation of Initial Miscible Time*  3.2.1. Solvent Accessible SURFACE Area (SASA) Analysis To confirm the initial miscible time discussed in previous parts, solvent accessible surface area (SASA) was calculated. SASA represents the hydrophobic, hydrophilic, and Based on similar treatments, crude oil atoms passing through the initial interface with time evolution (Figure S1) and its first-order variance (Figure S1b) can also be drawn. Figure 3 and Figure S1 show that the changes in CO<sup>2</sup> and crude oil were quite similar. Combined with the analyses above, we can preliminarily conclude that 4 ns is the initial miscible time, and is also the time when the CO2–oil system achieved miscibility.

#### total solvent accessible surface area for each component of the simulated system. Figure 4 *3.2. Reconfirmation of Initial Miscible Time*

*Molecules* **2021**, *26*, x FOR PEER REVIEW 5 of 14

space of the box along the Z-direction, as shown in Figure 2d.

#### shows the change in the hydrophilic area of CO2 from 0 ns to 10 ns and went through 3.2.1. Solvent Accessible SURFACE Area (SASA) Analysis

roughly three processes: (i) At the beginning, SASA increased rapidly and reached the highest point (from C1 to C2); (ii) SASA dropped to the lowest point in a short period of time (from C2 to C3); and (iii) SASA gradually rose to basic stability and attained a state of dynamic balance (from C3 to C4, and the time was after about 4 ns). To confirm the initial miscible time discussed in previous parts, solvent accessible surface area (SASA) was calculated. SASA represents the hydrophobic, hydrophilic, and total solvent accessible surface area for each component of the simulated system. Figure 4 shows the change in the hydrophilic area of CO<sup>2</sup> from 0 ns to 10 ns and went through roughly three processes: (i) At the beginning, SASA increased rapidly and reached the highest point (from C1 to C2); (ii) SASA dropped to the lowest point in a short period of time (from C2 to C3); and (iii) SASA gradually rose to basic stability and attained a state of dynamic balance (from C3 to C4, and the time was after about 4 ns).

> C1 and C2 are adjustments to the initial configuration in the molecular dynamics NPT ensemble, which was not our focus. This can be contributed to the molecular dynamic method being able to readjust the molecular conformation in the model under the NPT ensemble, and we focused more on the change in conformation after being readjusted. With

the blending of CO<sup>2</sup> and crude oil phases, both gradually achieved the best coexistence state (after C4).

Similarly, the changing trend of hydrophobic surface area of crude oil can be obtained by the same method. As shown in Figure 4b, a similar SASA change was observed in which the area increased and then decreased rapidly with time evolution (from O1 to O3 in Figure 4b). From Figure 4, it is reasonable to select 4 ns as the initial miscible time, and the data after 4 ns can be used to discuss the miscibility. *Molecules* **2021**, *26*, x FOR PEER REVIEW 6 of 14

**Figure 4.** SASA analysis for CO2 (**a**) and crude oil (**b**). **Figure 4.** SASA analysis for CO<sup>2</sup> (**a**) and crude oil (**b**).

### 3.2.2. Root Mean Square Deviation (RMSD) Analysis

best coexistence state (after C4).

C1 and C2 are adjustments to the initial configuration in the molecular dynamics NPT ensemble, which was not our focus. This can be contributed to the molecular dynamic method being able to readjust the molecular conformation in the model under Root mean square deviation (RMSD) compares each molecular structure in the simulation from the trajectory to the initial reference structure, reflects the change in its conformation, and is calculated by Equation (1).

$$\text{RMSD} = \sqrt{\frac{1}{N} \sum\_{i=1}^{N} \left( |r\_i(t) - r\_i(0)| \right)^2} \tag{1}$$

obtained by the same method. As shown in Figure 4b, a similar SASA change was observed in which the area increased and then decreased rapidly with time evolution (from O1 to O3 in Figure 4b). From Figure 4, it is reasonable to select 4 ns as the initial miscible time, and the data after 4 ns can be used to discuss the miscibility. 3.2.2. Root Mean Square Deviation (RMSD) Analysis where *N* is the total number of atoms (CO2/crude oil); and *r<sup>i</sup>* (0) and *r<sup>i</sup>* (t) are the initial position and the position of atom *i* at time *t.* Figure 5 displays the RMSD of CO2/crude oil during NPT ensemble as a function of time. It is interesting to note that CO<sup>2</sup> has a higher RMSD value than crude oil at the beginning, which indicates that CO<sup>2</sup> has better mobility. From 4 ns to 10 ns, the RMSD of CO2/crude oil in the box fluctuated with time evolution. Both were gathered around 4 nm of RMSD, which signifies that the system achieved equilibrium after 4 ns. *Molecules* **2021**, *26*, x FOR PEER REVIEW 7 of 14

> and represent the miscibility process between CO2 and crude oil phases. Interaction energy is a type of non-bonding interaction including long-range Coulomb interaction and short-range van der Waals interaction. As shown in Figure 6, the system was dominated by van der Waals interaction, while Coulomb interaction accounts for only about one-tenth of the former. This is because both CO2 and crude oil are non-polar molecules and do not have forces such as strong hydrogen bonding interaction. The intermolecular forces are mainly dispersive forces. The dispersion forces increased with time evolution, and the van der Waals potential energy and the total intermolecular

> When the system achieved equilibrium, the total interaction energy between CO2 and crude oil also reached its maximum and remained dynamically stable. Figure 6 clearly indicates that van der Waals interaction and Coulomb interaction both remained stable

> Once the initial miscible time system at 333 K and 10 MPa has been successfully determined, the number of CO2 atoms in the lower half of the box and crude oil atoms in

Both were gathered around 4 nm of RMSD, which signifies that the system achieved equilibrium after 4 ns. **Figure 5.** RMSD analysis for CO2 and crude oil. the upper half of the box after 4 ns were taken as the arithmetic mean respectively. It needs **Figure 5.** RMSD analysis for CO2 and crude oil.

potential energy increased accordingly.

**Figure 6.** Interaction energy analysis.

*3.3. Acquisition of MMP* 

after 4 ns.

3.2.3. Interaction Energy Analysis

### 3.2.3. Interaction Energy Analysis dominated by van der Waals interaction, while Coulomb interaction accounts for only

3.2.3. Interaction Energy Analysis

**Figure 5.** RMSD analysis for CO2 and crude oil.

*Molecules* **2021**, *26*, x FOR PEER REVIEW 7 of 14

The energy changes can reveal the changes in conformation in the simulated system and represent the miscibility process between CO<sup>2</sup> and crude oil phases. Interaction energy is a type of non-bonding interaction including long-range Coulomb interaction and shortrange van der Waals interaction. As shown in Figure 6, the system was dominated by van der Waals interaction, while Coulomb interaction accounts for only about one-tenth of the former. This is because both CO<sup>2</sup> and crude oil are non-polar molecules and do not have forces such as strong hydrogen bonding interaction. The intermolecular forces are mainly dispersive forces. The dispersion forces increased with time evolution, and the van der Waals potential energy and the total intermolecular potential energy increased accordingly. about one-tenth of the former. This is because both CO2 and crude oil are non-polar molecules and do not have forces such as strong hydrogen bonding interaction. The intermolecular forces are mainly dispersive forces. The dispersion forces increased with time evolution, and the van der Waals potential energy and the total intermolecular potential energy increased accordingly. When the system achieved equilibrium, the total interaction energy between CO2 and crude oil also reached its maximum and remained dynamically stable. Figure 6 clearly indicates that van der Waals interaction and Coulomb interaction both remained stable after 4 ns.

The energy changes can reveal the changes in conformation in the simulated system and represent the miscibility process between CO2 and crude oil phases. Interaction energy is a type of non-bonding interaction including long-range Coulomb interaction and short-range van der Waals interaction. As shown in Figure 6, the system was

**Figure 6.** Interaction energy analysis. **Figure 6.** Interaction energy analysis.

*3.3. Acquisition of MMP*  Once the initial miscible time system at 333 K and 10 MPa has been successfully determined, the number of CO2 atoms in the lower half of the box and crude oil atoms in When the system achieved equilibrium, the total interaction energy between CO<sup>2</sup> and crude oil also reached its maximum and remained dynamically stable. Figure 6 clearly indicates that van der Waals interaction and Coulomb interaction both remained stable after 4 ns.

### the upper half of the box after 4 ns were taken as the arithmetic mean respectively. It needs *3.3. Acquisition of MMP*

Once the initial miscible time system at 333 K and 10 MPa has been successfully determined, the number of CO<sup>2</sup> atoms in the lower half of the box and crude oil atoms in the upper half of the box after 4 ns were taken as the arithmetic mean respectively. It needs to point out that the number of CO<sup>2</sup> molecules under different pressures are different for 333 K system (Table S1). In order to reflect the general laws, the ratio of mean value to their respective total number of CO<sup>2</sup> and crude oil atoms in the box was calculated. Similarly, the ratio of CO2/crude oil passing through the initial interface to their respective totals at 333 K for 15 MPa, 20 MPa, 25 MPa, 30 MPa, and 35 MPa were also calculated, as shown in Table 1.

**Table 1.** The ratio of CO<sup>2</sup> and crude oil atoms passing through the initial interface to their respective totals.


From 10 to 35 MPa, the data of ratio decreased first and then became stable. We believe that this is because the system reached its peak pressure at 333 K. When the system exceeded this pressure, the additional simulation will not affect the value of each ratio. Therefore, the pressure is the theoretical MMP at 333 K. exceeded this pressure, the additional simulation will not affect the value of each ratio. Therefore, the pressure is the theoretical MMP at 333 K. For the sake of confirming the MMP, we handled the data according to its regularity. The first three decreasing points were fitted linearly, representing the systems before

From 10 to 35 MPa, the data of ratio decreased first and then became stable. We believe that this is because the system reached its peak pressure at 333 K. When the system

to point out that the number of CO2 molecules under different pressures are different for 333 K system (Table S1). In order to reflect the general laws, the ratio of mean value to their respective total number of CO2 and crude oil atoms in the box was calculated. Similarly, the ratio of CO2/crude oil passing through the initial interface to their respective totals at 333 K for 15 MPa, 20 MPa, 25 MPa, 30 MPa, and 35 MPa were also calculated, as

 **CO2 Crude Oil** 

10 MPa 0.0292 0.0300 15 MPa 0.0231 0.0236 20 MPa 0.0219 0.0221 25 MPa 0.0216 0.0217 30 MPa 0.0204 0.0207 35 MPa 0.0205 0.0208

**Table 1.** The ratio of CO2 and crude oil atoms passing through the initial interface to their

For the sake of confirming the MMP, we handled the data according to its regularity. The first three decreasing points were fitted linearly, representing the systems before MMP, and an equation in the form of y = kx + b was obtained. The last three nearly equal points were regarded as stable points, representing the systems after MMP, thus, another equation of y = x can be acquired by taking their arithmetic mean. We can subsequently obtain an intersection point as a consequence of simultaneous equations, and the abscissa corresponding to this point is the exact MMP at 333 K. As shown in Figure 7a,b, it was 20.31 MPa for CO<sup>2</sup> and 20.21 MPa for crude oil. MMP, and an equation in the form of y = kx + b was obtained. The last three nearly equal points were regarded as stable points, representing the systems after MMP, thus, another equation of y = x can be acquired by taking their arithmetic mean. We can subsequently obtain an intersection point as a consequence of simultaneous equations, and the abscissa corresponding to this point is the exact MMP at 333 K. As shown in Figure 7a,b, it was 20.31 MPa for CO2 and 20.21 MPa for crude oil.

*Molecules* **2021**, *26*, x FOR PEER REVIEW 8 of 14

shown in Table 1.

respective totals.

**Figure 7.** Acquisition of MMP (333 K) for CO2 (**a**) and crude oil (**b**). **Figure 7.** Acquisition of MMP (333 K) for CO<sup>2</sup> (**a**) and crude oil (**b**).

### *3.4. MMP in Different Temperature Systems*

*3.4. MMP in Different Temperature Systems*  We continued to simulate and analyze the data under the condition of 343 K, 353 K, 363 K, and 373 K at 10 MPa, 15 MPa, 20 MPa, 25 MPa, 30 MPa, and 35 MPa, respectively. It is worth pointing out that the density of CO2 varies greatly at different temperatures and pressures, therefore we computed the number of CO2 molecules under different conditions. The amount of CO2 molecules added to each simulation system are listed in Table S1. Afterward, a summary of initial miscible time in different systems can be We continued to simulate and analyze the data under the condition of 343 K, 353 K, 363 K, and 373 K at 10 MPa, 15 MPa, 20 MPa, 25 MPa, 30 MPa, and 35 MPa, respectively. It is worth pointing out that the density of CO<sup>2</sup> varies greatly at different temperatures and pressures, therefore we computed the number of CO<sup>2</sup> molecules under different conditions. The amount of CO<sup>2</sup> molecules added to each simulation system are listed in Table S1. Afterward, a summary of initial miscible time in different systems can be obtained according to the methods in Sections 3.1 and 3.2, as listed in Table 2. Table S2 summarized all the integral values in this study.


**Table 2.** Summary of initial miscible time (ns) in different systems.

The ratio of CO<sup>2</sup> and crude oil atom numbers that passed through the initial interface to their respective totals when they achieved miscibility can be obtained. Consequently, MMP of 343 K, 353 K, 363 K, and 373 K were obtained according to the method described in Section 3.4 by plotting and curve fitting (Figure S2). Table 3 summarizes the results.


**Table 3.** Summary of MMP (MPa) obtained from CO2/crude oil in different systems.

### *3.5. Model Assessment*

We fitted the MMP obtained from CO2/crude oil to TR, respectively, and obtained two prediction equations (Figure 8). It can be compared with the experimental results to check the predictive performance of the model. Recently, Yu et al. used a combination method of slim-tube experiments and interfacial tension (IFT) to perform MMP measurements on tight oil from the Long Dong region of the Ordos Basin. This method has higher credibility than slim-tube experiments [51]. Afterward, we compared our model with several famous empirical correlations to illustrate its accuracy by employing the experimental method proposed by Yu et al. as the benchmark. Table 4 reports the relative error. Details of these empirical correlations are summarized in Table S3. *Molecules* **2021**, *26*, x FOR PEER REVIEW 10 of 14

**Figure 8.** Relationships between TR and MMP for CO2 (**a**) and crude oil (**b**). **Figure 8.** Relationships between T<sup>R</sup> and MMP for CO<sup>2</sup> (**a**) and crude oil (**b**).


**Table 4.** Summary of MMP (MPa) and relative error predicted by experimental and different empirical correlations. **Table 4.** Summary of MMP (MPa) and relative error predicted by experimental and different empirical correlations.

> *3.6. Comparison of MMP Predicted by CO2 and Crude Oil*  The relationships between MMP predicted by CO2 and crude oil can be compared. It is more intuitive to reflect the data in Table 3 to Figure 9. In Figure 9, the blue line represents the curve whose analytical formula is y = x, and the red line is the fitting curve The relative error obtained from crude oil was similar to Lee [52], and the CO<sup>2</sup> relative error was similar to Alston et al. [53]. The equation proposed by Shokir [54] was based on an alternating conditional expectation algorithm, and had a relative error of 11.89 %. The model of Emera and Sarma [25] can be employed to calculate the MMP of impure CO<sup>2</sup> injection, but has poor accuracy. Beyond that, the performances of Cronquist [55],

> for the data. It can be found that the MMPs predicted by CO2 and crude oil were

values obtained from different reference materials (CO2/crude oil) were not identical as

there was a slight difference between them (i.e., an included angle of about 8°).

**Figure 9.** Comparison of MMP predicted by CO2 and crude oil.

Yellig and Metcalfe [57] 1 16.55 27.18

Glaso [56], and Yellig and Metcalfe [57] were also unsatisfactory. The overall results can prove that even if only the influencing factor of T<sup>R</sup> is considered, the model proposed in this study had satisfactory prediction accuracy. Glaso [56] 2 27.60 21.41 Yellig and Metcalfe [57] 1 16.55 27.18

### *3.6. Comparison of MMP Predicted by CO<sup>2</sup> and Crude Oil* The relationships between MMP predicted by CO2 and crude oil can be compared. It

*3.6. Comparison of MMP Predicted by CO2 and Crude Oil* 

*Molecules* **2021**, *26*, x FOR PEER REVIEW 10 of 14

**Figure 8.** Relationships between TR and MMP for CO2 (**a**) and crude oil (**b**).

**Model Number of Parameters Predicted MMP (MPa) Relative Error (%)** 

**Table 4.** Summary of MMP (MPa) and relative error predicted by experimental and different empirical correlations.

Yu et al. [51] - 22.75 - CO2 (this study) 1 19.63 13.71 Crude Oil (this study) 1 20.81 8.53

Lee [52] 1 20.84 8.32 Alston et al. [53] 4 19.72 13.22

Cronquist [55] 3 26.59 16.96

The relationships between MMP predicted by CO<sup>2</sup> and crude oil can be compared. It is more intuitive to reflect the data in Table 3 to Figure 9. In Figure 9, the blue line represents the curve whose analytical formula is y = x, and the red line is the fitting curve for the data. It can be found that the MMPs predicted by CO<sup>2</sup> and crude oil were approximately linear (R<sup>2</sup> = 0.955). Furthermore, in the same simulation system, the MMP values obtained from different reference materials (CO2/crude oil) were not identical as there was a slight difference between them (i.e., an included angle of about 8◦ ). is more intuitive to reflect the data in Table 3 to Figure 9. In Figure 9, the blue line represents the curve whose analytical formula is y = x, and the red line is the fitting curve for the data. It can be found that the MMPs predicted by CO2 and crude oil were approximately linear (R2 = 0.955). Furthermore, in the same simulation system, the MMP values obtained from different reference materials (CO2/crude oil) were not identical as there was a slight difference between them (i.e., an included angle of about 8°).

**Figure 9.** Comparison of MMP predicted by CO2 and crude oil. **Figure 9.** Comparison of MMP predicted by CO<sup>2</sup> and crude oil.
