Mechanisms for kerogen wettability transition from water-wet to CO2-wet: Implications for CO2 sequestration

Geological CO2 sequestration (GCS) is an essential building block of the global strategy to alleviate greenhouse gas emissions and mitigate the climate change. Injecting CO2 into the shale formations can not only reduce carbon emissions but also enhance oil recovery (EOR). Rock wettability is of great importance to CO2 storage as it determines the efficiency of structural and residual trapping of CO2 and plays a crucial role in CO2-EOR. In this work, molecular dynamics (MD) simulations are adopted to investigate the CO2-H2O-kerogen systems under various CO2 pressures. In a vacuum or under low CO2 pressures, kerogen surface is weakly water-wet thanks to the hydrogen bonding between H2O and kerogen. As CO2 pressure increases, kerogen wettability shifts from water-wet to CO2-wet, because more CO2 molecules accumulate at the H2O-kerogen interface and a distinct CO2 thin film emerges. Density functional theory (DFT) calculations reveal that the O-containing functional groups preferably adsorb H2O molecules over CO2 through hydrogen bonding, which is responsible for the weakly water-wet tendency at low CO2 pressures. In contrast, the carbon skeleton of kerogen exhibits a stronger affinity to CO2, leading to the formation of CO2 thin film on the kerogen surface. The CO2 crowding close to the kerogen surface at high CO2 pressures gives rise to the CO2-wet state. This study provides, for the first time, the fundamental mechanism for the kerogen wettability transition from water-wet to CO2-wet. The work also indicates that wettability of the mature kerogen is more likely to be CO2-wet during GCS, which is unfavorable for capillary trapping of CO2, but is favorable for CO2-EOR.


Introduction
Since the Industrial Revolution, the atmospheric CO 2 concentration has been climbing, reaching the level of ~420 ppm as of May 2021, which is ~50% above the pre-industrial level found in 1850 [1]. As a notorious greenhouse gas, excess amount of atmospheric CO 2 can cause global warming. Geological carbon sequestration (GCS) is a promising way to alleviate CO 2 emissions and forms one of the essential building blocks to reach the net-zero carbon emission target [2]. Shale formations have been considered as one of the ideal sites for GCS [3,4]. The geological conditions of shale formations are optimal for CO 2 storage [5] and the highly impermeable clay-bearing caprocks can also effectively seal the sequestrated CO 2 [6]. Furthermore, as GCS is often associated with a high cost, CO 2 enhanced oil recovery (EOR) in shale formations can greatly reduce the financial burden [2].
There are four main geological CO 2 storage mechanisms: structural trapping, where a caprock acts as a seal barrier preventing CO 2 from migrating to shallower zones [7]; residual or capillary trapping, where CO 2 is immobilized by capillary forces in pores and narrow throat in formation rocks [8]; solubility or dissolution trapping, where CO 2 dissolves in the formation water [9]; mineral trapping, where CO 2 reacts with rocks and formation water, forming carbonate minerals [10]. Among these mechanisms, the structural and residual trapping mechanisms are dependent on rock wettability, in which the CO 2 -H 2 O-rock contact angle plays a crucial role [11]. In the residual trapping, a large amount of free CO 2 gas is trapped by a high capillary force which is strongly dependent on the CO 2 -water-rock contact angle, while a waterwet rock is generally favored [12][13][14][15].
Shale media consist of inorganic and organic matters [16]. The inorganic matters include clay, quartz, and carbonate, and are generally considered hydrophilic [17]. The organic matters mainly consist of kerogen with its content up to 10 wt% [18]. Kerogen contains a considerable number of nano-scale pores with some pore throats size as small as 2 nm or less [19]. Such a narrow pore throat can serve as a prime sealing cap associated with an ultra-high capillary pressure for effective and stable structural and residual trapping. While the conventional wisdom is that kerogen is hydrophobic [20,21], a few recent experimental and simulation studies have shown that kerogen can be weakly water-wet [22][23][24][25], which might be attributed to the presence of heteroatoms (O, N, and S) in the kerogen structure [5]. Therefore, understanding about the CO 2 -water-kerogen contact angle plays an important role in GCS in shale formations. Besides, the CO 2 -waterkerogen contact angle is also closely related to Water alternating gas (WAG) flooding in shale formations [26] for EOR.
The CO 2 -water-rock contact angle is dependent on temperature [27][28][29], pressure [27,[30][31][32], rock type [33], pore structure and connectivity [34][35][36], and salt concentration [28,37,38], among other factors. Chiquet and coworkers [30] measured CO 2 -water-quartz and CO 2water-mica contact angles under different pressures, and found that as pressure increased from 0 to 11 MPa at 288 K, large increases in the contact angles were observed in the order of 15-25 • for quartz and 40-50 • for mica. As a result, these two substrates shifted from a waterwet state under a low CO 2 pressure to an intermediate-wet state under a high CO 2 pressure at 288 K. They suggested that the wettability alteration was caused by the pH reduction as well as the increase in CO 2 pressure. Similar phenomena of rock wettability transition were also reported by other experimental studies in the CO 2 -water-quartz [31,39] and the CO 2 -water-mica [27] systems. Much less attention has been paid to the contact angles of CO 2 -water-organic rock systems. Arif et al. [40] measured the contact angles of CO 2 -water-coal systems and reported that when the pressure increased from 0.1 MPa to 20 MPa at 323 K, the advancing contact angle increased from 47 • to 141 • and the receding contact angle increased from 42 • to 129 • . They explained that the increased contact angle was related to the increased CO 2 adsorption under higher pressure. As for the CO 2 -water-kerogen system, Ho et al. [25] reported that the water contact angle on kerogen surface was 42.8 • in the absence of CO 2 , while the contact angle was 180 • when the CO 2 pressure was 200 atm at 300 K, indicating that the kerogen wettability shifted from water-wet to CO 2 -wet. They observed that CO 2 formed a thin film on the kerogen surface which might facilitate the wettability change. However, their work focused on CO 2 effects on the water flow in the kerogen nanochannel and paid little attention to the CO 2 -induced wettability transition mechanism. While these works reported contact angles in CO 2 -water-rock systems over a wide range of pressure, temperature and rock chemistry, there is a lack of direct investigations into the interactions between water and kerogen in the absence/presence of CO 2 , and the underlying mechanisms of kerogen wettability transition remain elusive.
A comprehensive understanding of the wettability in CO 2 -waterkerogen systems is needed for geological CO 2 sequestration and CO 2 -EOR as kerogen is the main constituent of organic matters in shale. In this work, we study interfacial phenomena and evaluate contact angles in CO 2 -water-kerogen systems using molecular dynamics (MD) simulations under typical GCS conditions. The contact angle in the CO 2 -H 2 Okerogen system increases until reaching 180 • as CO 2 pressure increases. The interfacial configurations, density profiles and hydrogen bond densities are carefully analyzed. Moreover, the effect of kerogen heteroatoms on the CO 2 -induced contact angle transition is explained by the first-principle density functional theory (DFT) calculations. The work provides a fundamental understanding about CO 2 -H 2 O-kerogen interactions and reveals the underlying mechanism of wettability variability from sub-atomic and molecular perspectives.

Molecular dynamics simulations
The schematic representation of the CO 2 -H 2 O-kerogen system is shown in Fig. 1. A typical simulation system consists of a cylindrical water droplet residing on the kerogen surface in a CO 2 environment. The simulation box has the dimensions of 110 × 135 × 172 Å 3 . The simulation domain is tested to be large enough to eliminate any finite size effect. The cylindrical droplet geometry is employed to avoid the line tension effect [41]. The kerogen surface slab is constructed by MD simulations following the procedures reported in our previous works [42][43][44][45]. The Type II-D kerogen molecular model proposed by Ungerer et al. [46] is chosen to represent mature kerogens in shale formations. During the simulation, the kerogen structures are modelled using the consistent valence force field (CVFF) [47]. The final size of the kerogen slab is 110 × 135 × 36.5 Å 3 . The density of the simulated kerogen is 1.19 g/cm 3 , which falls within the range of experimental values of type II kerogen of 1.18-1.35 g/cm 3 [48]. Then, a water droplet consisting of 8,000 H 2 O molecules is placed on top of the kerogen surface. Fig. S1 in the Supporting Information (SI) shows the initial configuration of the H 2 O-kerogen system in the absence/presence of CO 2 . The number of CO 2 molecules in the surrounding CO 2 phase ranges from 0 to 23,300 to represent various pressure conditions ranging from 0 to 60 MPa. The CO 2 density far away from the H 2 O droplet and kerogen surface is used to dictate the pressure P by comparing to the NIST Chemistry Webbook [49]. Table 1 summarizes the number of H 2 O and CO 2 molecules in the simulation systems. For each case, the values of pressure and the uncertainties were the average value and the standard deviation, respectively, calculated from 5 independent estimates obtained by dividing the 5 ns production run into segments of 1 ns.
Throughout the simulations, the kerogen surface is held stationary, while H 2 O and CO 2 molecules can move freely. The H 2 O molecules are described by the SPC/E model [50], while we use the rigid TraPPE model to represent the CO 2 molecules [51]. The combination of SPC/E-TraPPE force fields has been proven to be very accurate for describing CO 2 -H 2 O interaction in the temperature range 323.15-478.15 K and pressures up to 100 MPa [52], in comparison with the experimental data [53]. The Lorentz-Berthelot mixing rules [54] are employed to calculate the interactions between unlike atoms. The non-bonded interactions are represented by the pairwise Coulomb potential and Lennard-Jones (LJ) potential. A cutoff radius of 1 nm is used to compute the LJ interactions with the analytical tail corrections [55]. The long-range electrostatic potential is computed using the PPPM integrations [56]. Periodic boundary conditions are applied in both the x-and y-directions. In the zdirection, a virtual wall is placed at the top of the simulation box which is far away from the H 2 O droplet and kerogen surface to confine the CO 2 molecules. The fluid-virtual wall interaction is given by [57,58], where z denotes the distance between the fluid molecule and the wall, ε sf and σ sf are obtained from the Lorentz-Berthelot mixing rules of energy and size parameters of CO 2 molecules, respectively. The potential was cut and shifted at its minimum, to ensure that the virtual wall only exerts repulsive force to CO 2 molecules. A vacuum with the same sizes in x-and y-directions and three times of the length of the simulation box in the zdirection is inserted to avoid the artificial influence from the periodic images in the z-direction [55,59]. The entire simulation box including the vacuum space is shown in Fig. S2. All MD simulations are carried out using the open-source package LAMMPS [60]. The time step is 1 fs. All systems are initially equilibrated for 10 ns in the NVT ensemble and the temperature is maintained at 330 K by using the Nosé-Hoover thermostat [61] with a damping constant of 100 fs. The system is considered to reach equilibrium when little fluctuation can be observed in the density distributions. Then, production runs of 5 ns are conducted and the trajectory snapshots are analyzed to obtain the CO 2 -H 2 O-kerogen contact angle. Further details regarding the contact angle measurement can be found in the SI.

DFT calculations
We use DFT to investigate the interactions among CO 2 , H 2 O and the surface functional groups on kerogen. Competitive adsorption of H 2 O and CO 2 on different functional groups are simulated and the role of the main functional groups in the H 2 O-CO 2 co-adsorption is studied. The first-principle DFT with dispersion correction of D3 scheme (DFT-D) calculations is implemented in Vienna ab-initio Simulation Package (VASP) [62,63]. All the geometries for modelling are established in MedeA 3.1 from the Materials Design. The projector augmented-wave (PAW) [64] method and the generalized-gradient-approximation (GGA) [65] exchange-correlation functional in the form of Perdew-Burke-Ernzerhof (PBE) [66] are employed in all the calculations. The electronic wave function is expanded on a plane-wave basis with the energy cutoff of 600 eV to converge the relevant quantities. The Brillouin zone is sampled using a 2 × 2 × 1 Monkhorst-Pack k-point with a smearing of 0.1 eV. The self-consistent field (SCF) tolerance is set as 10 − 6 eV/atom. The entire calculation is performed with a convergence threshold of 0.03 eV/Å on the maximum force. No symmetry constraint is used for any modelling. Graphene-based models derived from optimized graphite cell (see Fig. S3) are adopted to simulate the kerogen facets; a 4 × 4 supercell model is used for the graphene and a 6 × 6 supercell model is used for the O-doped graphene, in which the dangling bonds are passivated by H. A 15 Å vacuum region is created above the top of the facets. The adsorption energy is calculated as [67], where E total is the total energy of the substrate with the adsorbed gas molecule, E substrate denotes the energy of the bare substrate and E adsorbate is the energy of the isolated gas molecule. It is widely accepted that the most negative adsorption energy means the most stable configuration. The energy defined by the following equation is used to evaluate the total interactions of the adsorbates, with the assumption that the variation of adsorption energies over different adsorption sites is negligible.
where E tot is the total co-adsorption energy for the whole system, E CO2 and E H2O are the total adsorption energies for CO 2 and H 2 O, respectively.

Results and discussion
In this section, we first calculate the CO 2 -H 2 O-kerogen contact angle in the absence or presence of CO 2 . The effect of CO 2 pressure on the contact angles is analyzed in detail. Then, the underlying mechanisms of kerogen wettability transition induced by the CO 2 pressure increase are discussed via MD and DFT simulations.

Contact angles on kerogen surface
Two examples of representative equilibrium snapshots and the corresponding two-dimensional (2-D) density contours of water droplets at 0 MPa and 11.24 MPa and 330 K are plotted in Fig. 2a and b. The 2-D density contours of H 2 O and CO 2 under the other pressures are given in Fig. S4. When P⩽4.66 MPa, the water droplet resides on the kerogen surface and two distinct H 2 O adsorption layers are observed on the kerogen surface. The formation of H 2 O layering structures can be attributed to the hydrogen bonding between H 2 O and kerogen heteroatoms as we will discuss below as well as the entropic effect [68]. As P increases, the two layers become less prominent and the H 2 O density distributions close to the kerogen surface become more uniform. At P = 44.40 MPa, the water droplet is completely detached from the kerogen surface. Fig. S4d and f clearly show uniform CO 2 distributions away from the water droplet with its density close to the bulk one. However, in the vicinity of the kerogen surface, significant CO 2 adsorption layers can be observed on both sides of the water droplet. As P increases, the CO 2 density at the H 2 O-kerogen interface gradually increases and a distinct CO 2 film emerges at 44.40 MPa, which separates the water droplet from the kerogen surface.
To illustrate the dependence of CO 2 -H 2 O-kerogen contact angle θ on P, θ versus CO 2 pressure is plotted in Fig. 2c. The kerogen wettability can   [69]. Combined with the analysis on CO 2 density contour plots ( Fig. S4d and f), the effect of CO 2 pressure on θ can be attributed to the accumulation of CO 2 molecules at the H 2 O-kerogen interfaces.  Fig. 3a, we depict the H 2 O and CO 2 density profiles along the center of the water droplet at 4.66 MPa. The density profiles are averaged over a 12 Å-wide strip in the center along the x-direction as illustrated in Fig. 3. The topmost atom of the kerogen slab is defined as the origin in the z-direction. The uneven kerogen surface is responsible for the non-zero density at z < 0. Close to the kerogen surface, two prominent H 2 O density peaks can be observed at z = 0.24 nm and z = 0.54 nm with their values at 1.15 and 1.03 g/cm 3 , respectively, in line with the layering structures observed in Fig. 2. The H 2 O density regresses to its bulk value (0.97 g/cm 3 ) inside the droplet. On the other hand, CO 2 density profiles render two distinct peaks: one located at the H 2 O-kerogen interface; the other one at the H 2 O-CO 2 interface. The peak in the vicinity to the kerogen surface indicates the accumulation of CO 2 due to the CO 2 -kerogen interactions which are largely hydrophobic interactions. In addition, CO 2 is dissolved within the H 2 O droplet. Away from the H 2 O droplet, CO 2 density regresses to its bulk value. To better understand the structural properties of the CO 2 -H 2 O-kerogen system, in Fig. 3b and Fig. 3c, we present the corresponding snapshots showing the entire system and a close-up view of the H 2 O-kerogen interface region. They show that a number of CO 2 molecules accumulate at the H 2 Okerogen interface. Fig. 4(a) Mass density profiles of H 2 O and CO 2 averaged over the 12 Å-wide strip normal to the kerogen surface at 15.95 MPa. The snapshot of (b) the entire system; (c) a close-up view of H 2 O-kerogen interface region. The color code is the same as Fig. 1.

Density distributions
The thin CO 2 film at the H 2 O-kerogen interface continues to grow as CO 2 pressure reaches 34.26 MPa as shown in Fig. 5a. Comparing to Fig. 4a, the peak value in CO 2 density profiles in the vicinity of kerogen surface further increases. As a result, H 2 O is depleted from the kerogen surface. The enrichment of CO 2 at the CO 2 -H 2 O interface disappears at high pressures, in line with previous work [70]. As shown in Fig. 5b, the water droplet is almost completely detached from the kerogen surface and a distinct CO 2 film can be observed. In contrast to Fig. 3c and Fig. 4c, more CO 2 molecules accumulate at the H 2 O-kerogen interface.
To reveal the pressure effect on the shape of water droplet, the evolution of droplet contour profiles, defined as the boundary between  the droplet and the surrounding where the H 2 O mass density equals the CO 2 mass density, is depicted in Fig. 6. It is apparent that as the CO 2 pressure increases, the radius of the water droplet decreases. During this process, the kerogen surface turns from H 2 O-wet to CO 2 -wet. To quantify the relation between the CO 2 pressure and the CO 2 thin film at the H 2 O-kerogen interface, the CO 2 film thickness δ is computed. Due to the roughness of the kerogen surface, it is difficult to distinguish the actual δ. Therefore, we define the intersection of the contour profile and the dashed line in Fig. 6 at 0 MPa as the zero reference (δ = 0). The distance from the intersection to δ = 0 is regarded as an indicator of δ. δ as a function of the CO 2 pressure is displayed in the inset of Fig. 6. The CO 2 film thickness first rapidly increases to 0.30 nm at 0.65 MPa, then the increasing rate slightly declines until 15.95 MPa and the CO 2 film thickness eventually reaches a plateau when the CO 2 pressure further increases to 44.42 MPa. Increase in CO 2 film thickness can be attributed to increasing CO 2 adsorption on kerogen surface with increasing CO 2 pressure. The plateau occurs when a complete adsorption layer is formed. As the water droplet leaves the kerogen surface at a higher pressure indicating a completely CO 2 -wet state, we only report the film thickness until 44.42 MPa.  Fig. 4a, a sharp peak at z = 0.14 nm is observed in the CO 2 density profile, which highlights the strong accumulation of CO 2 at H 2 O-kerogen interfaces, consistent with the finding of Ho et al. [25]. As a result, the H 2 O layering structures in the vicinity of kerogen surface become less prominent. This phenomenon can be explained by the snapshots shown in Fig. 4b and c. On the other hand, the enrichment of CO 2 at CO 2 -H 2 O interface becomes less significant.  Fig. 1.   Fig. 6. Water droplet contour profiles as a function of CO 2 pressure. Inset: CO 2 film thickness as a function of CO 2 pressure. The dashed line represents the center of the water droplet.

Hydrogen bond analysis
To further reveal the changes in the water affinity to the kerogen surface with increasing CO 2 pressure, the hydrogen bonds between water and kerogen surface are analyzed. The geometric criterion proposed by Luzar and Chandler [71] is adopted. A hydrogen bond exists if the distance between the donor and acceptor is less than 0.35 nm with the hydrogen-donor-acceptor angle less than 30 • . A recent research [72] reported that the strongest interaction between H 2 O and kerogen is contributed by hydrogen bonds in which H 2 O serves as a donor and the heteroatoms (N, O, and S) in kerogen serve as an acceptor. The N atom shows the highest affinity to H 2 O, followed by O and S. In this work, the number of hydrogen bonds between H 2 O and the kerogen heteroatoms is calculated and averaged over the 5-ns production runs. As shown in Fig. 7, it is evident that as CO 2 pressure increases, the total number of hydrogen bonds drops gradually. The results demonstrate that the intrusion of CO 2 into the H 2 O-kerogen interface greatly diminishes the H 2 O-kerogen affinity. In all cases, the O atoms provide the largest fraction of hydrogen bonds (at least 77%), while the S atoms make little contribution. The large hydrogen bond number for O atoms may result from their abundance in kerogen as highlighted in the inset in Fig. 7. The number of O atoms is 2.25 times that of N atoms. The important role of O content in kerogen wettability was also reported by Hu et al. [21]. They found that kerogen wettability may shift from hydrocarbon-wet to water-wet by increasing the number of carbonyl (-C--O) groups on the graphene surface. With the O/C ratio at 20%, the graphene surface with heterogeneously-distributed carbonyl groups is completely water-wet. Herein, graphene and O-doped graphene models are chosen to simulate the representative carbon skeleton and O-containing functional groups on the kerogen surface, respectively. The graphene model is derived from a fully optimized graphite cell and the calculated lattice constants has been validated in Fig. S3. Firstly, we compare the adsorption energy of CO 2 and H 2 O molecules on the graphene surface (carbon skeleton of kerogen). The adsorption energy of H 2 O on the graphene surface is − 0.12 eV, which is consistent with the range of − 0.07 eV to − 0.13 eV found in previous studies [73][74][75][76][77]. The adsorption energy of CO 2 on graphene surface is − 0.15 eV, which is lower than that of H 2 O, suggesting that the adsorption of CO 2 on graphene surface is more favorable than H 2 O. Then, we investigate the interaction between H 2 O and graphene surface (carbon skeleton) in the presence of CO 2 . Fig. 8 shows the co-adsorption configurations for H 2 O and CO 2 on the graphene surface, and the interaction energy E int as defined in Eq. (3) is listed for each configuration. The E int value of one H 2 O molecule on the graphene surface is set at zero as a reference. The E int values for co-adsorption of an H 2 O-CO 2 mixture (− 0.11 eV) and H 2 O-2CO 2 mixture (− 0.24 eV) on the graphene surface are negative, arising from the interactions among CO 2 and H 2 O in both systems. In particular, new hydrogen bonds form between the CO 2 and H 2 O molecules. In addition, the presence of CO 2 molecules shifts the H 2 O molecule away from the graphene surface. The results imply that H 2 O adsorption on the graphene surface is weakened by the presence of CO 2 . For further verification, MD simulations are carried out to obtain the contact angles of water droplets on graphene surface under different CO 2 pressures as shown in Fig. S5. Under the same CO 2 pressure, the contact angle on the graphene surface is larger than that on the kerogen surface. The water droplet is completely detached from the graphene surface at 16.32 MPa, while the same phenomenon occurs at 44.42 MPa for the CO 2 -H 2 O-kerogen system. The results reveal that the aromatic carbon skeleton of kerogen preferably adsorbs CO 2 compared to H 2 O, and CO 2 can form a thin film to deplete H 2 O molecules from the kerogen surface, in line with our MD simulation results in Section 3.2.

Effect of surface functional groups
As discussed in Section 3.  Fig. 9. In the first scenario, two H 2 O molecules are located near the doped Oe site, leading to the formation of two hydrogen bonds: one is between the two H 2 O molecules; the other is between the H 2 O molecule and the doped Oe on the substrate, as shown in Fig. 9a. To illustrate the CO 2 effect, the adsorption of H 2 O on the O-doped graphene surface is modelled in the presence of CO 2 molecules, as shown in Fig. 9b and c. The original adsorption positions of the two H 2 O molecules in the absence of CO 2 are also displayed in green for comparison. As shown in Fig. 9b,   While this result is seemingly in a stark contrast to the role of CO 2 pressure on CO 2 -H 2 O-kerogen contact angle as discussed above, we note that the DFT simulations are conducted at the conditions in which there are a very small number of CO 2 molecules. It rather endorses the formation of hydrogen bonding between H 2 O and kerogen surface, and explains why kerogen still behaves as a weakly water-wet [22][23][24][25] or intermediate-wet substrate at a relatively low CO 2 pressure as shown in Fig. 2. On the other hand, as CO 2 pressure further increases, due to the strong interactions between CO 2 and carbon skeleton as shown in Fig. 8, more CO 2 molecules accumulate in the vicinity of kerogen surface. The formation and constant-growth of CO 2 thin film would influence the formation of hydrogen bonding between H 2 O and kerogen surface. This phenomenon can be attributed to the correlation effect among CO 2 , H 2 O and kerogen which can be fully captured by MD simulations. On the other hand, DFT simulations cannot grasp such effect due to high computational costs. The correlation effect stems from intermolecular interactions and molecular configurations and is important when a large number of fluid molecules gather at the interface [78]. For example, a specific fluid molecule distribution close to the substrate is dependent not only on fluid-surface interactions but also molecular configurations of the surrounding fluid. In other words, in a vacuum or at a low CO 2 pressure, the CO 2 Fig. S5. Under the same CO 2 pressure, the water contact angle on the graphene surface is larger than that on the O-doped graphene surface, and the H 2 O droplet is totally detached from the graphene surface at 16.32 MPa, while the same phenomenon occurs at 18.77 MPa for the Odoped graphene system. Besides, as shown in Fig. S5c, the density of the CO 2 adsorption layer near the graphene surface is larger than that of the O-doped graphene surface under 8.85 MPa. These results further verify the CO 2 pressure effects and reveal that CO 2 can potentially form a thin film with increasing pressure and repel H 2 O molecules from the kerogen surface.

Conclusions
In this work, the effect of CO 2 pressure on the CO 2 -H 2 O-kerogen contact angle is investigated through MD and DFT simulations. The simulation results show that the contact angle increases from 60.4 • to 180 • when the CO 2 pressure increases from 0 to 44.42 MPa. The presence of CO 2 results in a gradual transition in the kerogen wettability from H 2 O-wet to CO 2 -wet. At low CO 2 pressures, CO 2 and H 2 O co-adsorb at the solid-fluid interface. As CO 2 pressure increases, a distinct CO 2 film is observed at the H 2 O-kerogen interface, leading to a sharp decline in the number of hydrogen bonds between H 2 O and kerogen. The H 2 Okerogen interaction is dominated by the hydrogen bonding, while the Ocontaining functional groups on the kerogen surface contribute the majority of hydrogen bonds. According to the DFT calculations, the presence of CO 2 pushes the H 2 O molecule closer to the Oe group on the O-doped graphene substrate. However, on the graphene surface, CO 2 has a lower adsorption energy than H 2 O, while CO 2 can push H 2 O away from the substrate. Therefore, in a vacuum or at a low CO 2 pressure, the CO 2 -H 2 O-kerogen contact angle is mainly determined by H 2 O-kerogen and H 2 O-H 2 O interactions, while it is greatly affected by CO 2 crowding driven by the hydrophobic interactions at a high CO 2 pressure. Our work indicates that the highly-matured kerogen with a relatively low O content is prone to wettability change in the presence of CO 2 . The CO 2 -wet state at high pressures is unfavorable for CO 2 residual trapping, but is favorable for the water alternating gas process for EOR.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.