The pore model data for kaolinite used in this study were obtained from experiments (Bish and Dreele, 1989). The initial crystal cell belongs to the P1 space group, indicating a monoclinic crystal system. The molecular formula is Al4Si4O108, comprising a total of 34 atoms. The structure is classified as a layer-type mineral with TO-type octahedral layers. There is only one unit layer in the middle, and it does not contain any water ions or sodium ions. The lattice parameters are as follows: a = 0.514 nm, b = 0.893 nm, c = 0.737 nm, α = 91.8°, β = 104.5°, γ = 90.0°. The initial crystal cell used in the simulation is from the kaolinite crystal cell in the Materials Studio software database. After structural optimization, the relative errors between the crystal cell and experimental data are shown in Table 1. All relative errors are less than 0.5%. Therefore, it is evident that the chosen kaolinite pore structure model for simulation is reasonably accurate. The calculation of relative error is shown by Formula 1: δ R = X e − X s X e × 100 % (1)

Where, δ R is relative error, X e is experimental unit-cell, X s is simulate unit-cell.

The pore structure of kaolinite used in this study was constructed by expanding and combining crystal cells from the Materials Studio database. There are two types of pore surfaces, one is the (0 0 1) surface with hydrogen atoms at the end of the crystal cell, the other is the (0 0–1) surface without hydrogen atoms. The surface was extended using a supercell to create a surface with dimensions of 5 a × 4 b × 2 c, where the size of the XOY plane is 2.574 nm × 3.574 nm. This expansion aimed to increase the interaction positions between carbon dioxide and the kaolinite pore. Subsequently, layer-like pores were constructed to simulate micro/nano pores, There are three types of configurations for the inner surface of layered narrow pores, hydrogen containing and hydrogen containing surfaces, hydrogen free and hydrogen free surfaces, hydrogen containing and hydrogen free surfaces. In this article, the most common (0 0 1) hydrogen-containing surface is selected as the inner surface of the layered structure of kaolinite (Bin et al., 2018), two kaolinite layers were placed face-to-face parallel to the XOY plane, and different pore widths were achieved by adjusting the thickness of the vacuum layer, as illustrated in Figure 2.

The carbon dioxide model used in the simulation adopted the linear EPM2-type CO₂ model (Harris and Yung, 1995), widely applied in micro-scale simulations of carbon dioxide adsorption (Liao et al., 2022). The bond length of carbon dioxide is 0.149 nm, and the bond angle is 180°, representing a double-bond structure between carbon and oxygen.

The kaolinite pore model employs the Clayff force field, which has been demonstrated to be widely applicable for describing the interactions of clay minerals (Cygan et al., 2004; Cygan et al., 2009). The potential energy parameters for carbon dioxide adopt the EPM2 model parameters (Harris and Yung, 1995; Qin et al., 2008). The charges and Lennard-Jones potential energy parameters for kaolinite and carbon dioxide are presented in Table 2. The intra-molecular bond energy, bond angle energy, and dihedral angle energy are all described using harmonic functions. The non-bonded potential energy between molecules is described by the following Formula 2: U r i j = U L J + U C = 4 ε i j σ i j r i j 12 − σ i j r i j 6 + q i q j 4 π ε 0 r i j (2)

Where r i j is the distance between atoms i and j, U r i j is the non-bonded potential, U L J is the Lennard-Jones potential, U C   is the Coulombic interaction potential, q i and q j is the charge on atoms i and j, ε 0 is the permittivity of vacuum, ε i j is the Lennard-Jones potential well depth, and σ i j is the Lennard-Jones zero potential distance. The parameters ε i j and σ i j are determined using the Lorentz–Berthelot combining rules (Formula 3; Formula 4): σ i j = σ i + σ j 2 (3) ε i j = ε i ε j (4)

All molecular dynamics simulations in this study were conducted using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) software (Plimpton, 1995; Zhou et al., 2020; Liao et al., 2021; Meng et al., 2021). The simulations employed three-dimensional periodic boundary conditions, the long-range Coulomb interaction were treated using the PPPM summation (Brown et al., 2012), the calculated relative error is 10–4, utilizing the NVT ensemble, and the temperature control was achieved through the Nosé-Hoover algorithm (Evans et al., 1985). The intermolecular cutoff radius was set to 1.0 nm based on the findings by Liao et al., 2022, indicating that this cutoff radius is sufficient for calculating the non-bonded interactions of carbon dioxide molecules. Throughout the simulation process, the kaolinite pore was kept fixed to maintain stability, treated as a rigid body where the positions of reservoir pores do not update with increasing time steps. Consequently, the shape and size of the pore remain constant during the simulation.

After obtaining models of kaolinite pores with different widths and an ideal carbon dioxide molecule, the microscopic flow behavior of carbon dioxide in the nanoscale pores was constrained using the relevant force field parameters. Initially, carbon dioxide molecules of equal volume density were inserted into the vacuum region of the constructed pore, as illustrated in Figure 2. As the pressure of the carbon dioxide fluid is a function of pore width, the number of molecules, and temperature, varying the temperature results in different pressures for different numbers of carbon dioxide molecules in pores of different widths. Subsequently, a 6 ns equilibrium molecular dynamics (EMD) molecular dynamics simulation was performed on the entire constructed model, with a time step of 1 fs. This process achieves thermodynamic equilibrium of the system and allows carbon dioxide molecules in the micro/nano pores to converge to equilibrium pressure values. The number of inserted carbon dioxide molecules and the corresponding equilibrium pressure values after the system stabilizes are presented in Table 3. The last 1.5 ns of the 6 ns EMD simulation were selected for thermodynamic data collection of the pore fluid. Using the EMD method, the state and distribution characteristics of carbon dioxide within the pore can be obtained.

Upon achieving thermodynamic equilibrium using the EMD method, an acceleration parallel to the XOY surface of the pore was applied to each atom of the carbon dioxide fluid within the pore. To ensure that the molecular velocity overcomes thermal velocity in the specified direction, the simulation time for non-equilibrium molecular dynamics (NEMD) must be significantly longer than that for EMD. Therefore, a total of 12 ns of NEMD molecular dynamics simulation was conducted, with a time step of 1 fs. After the system reached dynamic equilibrium, data were collected every 100 steps, and the last 2 ns of the NEMD simulation were chosen for thermodynamic data collection and subsequent analysis of statistical results and trajectory files.

This study employs non-equilibrium molecular dynamics (NEMD) simulations to investigate the flow behavior of carbon dioxide in kaolinite pores under the application of a pressure gradient. Acceleration is applied to each atom of the carbon dioxide fluid to simulate the flow of carbon dioxide within the kaolinite pores induced by the pressure difference (Abbasi and Karimian, 2016). The acceleration applied by all carbon dioxide molecules in the simulation is 0.02 nm/ps2, which means the pressure difference is 42.7 MPa, The relationship between acceleration (a) and pressure difference (ΔP) is described by Formula 5: ∆ P = N m a y A x z (5)

Here, N is the number of particles in the pore, m is the particle mass, a y is the acceleration applied in the y-direction, and A x z is the cross-sectional area of the pore perpendicular to the acceleration in the y-direction.

Due to the well-developed interlayer pores in clay minerals, the surfaces of these micro/nano pores serve as the primary locations for fluid adsorption. According to the International Union of Pure and Applied Chemistry (IUPAC) definition (Sing et al., 1985), micro/nano pores are generally classified into three categories: micropores (less than 2 nm), mesopores (2 nm–50 nm), and macropores (greater than 50 nm). Among these, mesopores are most widely distributed in kaolinite. Therefore, this study focuses on mesoporous pores as the main research subject. Considering that with increasing mesopore width, the influence of the pore wall on the fluid inside the pore gradually weakens, and carbon dioxide within the pores is primarily in the free state, pore widths of 2 nm, 3 nm, and 5 nm are selected as the study objects. Assuming surface temperature is 6.85°C, surface pressure is 0.101 MPa, average geothermal gradient is 25°C/km, and hydrostatic pressure gradient is 0.0098 MPa/m (Zhang et al., 2015), different temperature and pressure conditions corresponding to burial depths of 1200 m, 2,400 m, 3,600 m, and 4,800 m were selected. These depths correspond to temperatures and pressures as follows: 1200 m: 310 K, 11.9 MPa; 2,400 m: 340 K, 23.6 MPa; 3,600 m: 370 K, 35.4 MPa; 4,800 m: 400 K, 47.1 MPa. Comparing the pressure-temperature distribution curve in Figure 3 with the assumed ideal conditions at different burial depths, a strong agreement in basic trends and patterns is observed, indicating the reliability of the simulation results. Therefore, this study aims to investigate the state of carbon dioxide within kaolinite pores and the impact of flow behavior induced by pressure gradients on the carbon dioxide fluid’s state at burial depths ranging from 1,200 m to 4,800 m.

This study employs the Equilibrium Molecular Dynamics (EMD) method to investigate the distribution of carbon dioxide (CO2) fluid within pores of 2 nm, 3 nm, and 5 nm diameters. The density distribution at a burial depth of 2.4 km is calculated. Figure 4A illustrates the schematic distribution of CO2 within the pores, while Figure 4B presents a density distribution contour map for a 3 nm pore width. Both axes are in Å, and the color scale is in g/cm3, representing different density distributions within the pores.

From Figure 3A, it is observed that the density distribution of CO₂ is similar for various pore widths. Within the kaolinite pores, a symmetrical distribution is evident, with two distinct adsorption peaks near the surface. This pattern aligns with previous studies by Zhang et al. (Zhang and Jin, 2022) on CO₂ adsorption on kaolinite surfaces. The first adsorption peak is attributed to van der Waals forces and hydrogen bonding between CO₂ molecules and hydroxyl groups on the kaolinite pore surface. The second peak primarily arises from van der Waals forces between CO₂ molecules, as a stable adsorption layer has already formed at the surface. Due to structural reasons, the van der Waals forces between nonpolar CO₂ molecules are significantly weaker than those between CO₂ and the pore surface, explaining the weaker intensity of the second adsorption peak compared to the first. The central region of the pore exhibits a relatively uniform density distribution, with CO₂ molecules mainly in a free state. Comparing the thickness of adsorption layers, the first layer’s thickness remains consistent, while the second layer’s thickness follows the order: 5 nm = 3 nm > 2 nm. This trend is attributed to the increased influence of the pore’s upper and lower surfaces on CO₂ molecules as the pore width decreases, making the formation of a saturated second adsorption layer more challenging.

Figure 3B illustrates that near the kaolinite surface, colors are predominantly red, yellow, and green, indicating higher CO₂ density compared to the central region. This is due to the weak polarity of CO₂ molecules, leading to van der Waals forces with polar hydroxyl groups on the kaolinite surface, and hydrogen bonding between oxygen atoms in CO₂ and hydrogen atoms in hydroxyl groups. In contrast, regions farther from the kaolinite surface appear in shades of blue, indicating lower overall CO₂ density. This is attributed to the diminishing influence of the kaolinite surface on CO₂ as distance increases, allowing more CO₂ molecules to diffuse freely in the central region.

To further investigate the density distribution of CO₂ at different burial depths, simulations were conducted for varying depths within kaolinite pores. Figure 5 compares the density distribution of CO₂ at different depths, revealing a decreasing trend in both adsorption peaks with increasing burial depth. The first adsorption peak experiences a more significant reduction, decreasing from 1.904 g/cm³ to 1.445 g/cm³, decrease of up to 24.1 percent. As burial depth increases, temperature rises, intensifying the thermal motion of CO₂ molecules. This results in more CO₂ molecules desorbing from the adsorption layer and diffusing towards the central region. The second adsorption layer and free-state CO₂ molecules, experiencing weaker confinement, frequently collide with the stable first adsorption layer, leading to a greater reduction in its density. The desorbed CO₂ molecules contribute to an increasing density trend in the central fluid of the pore, that’s an increase of almost 9 percent.

To further understand the factors influencing the carbon dioxide (CO₂) storage state with varying burial depths, this study employs the strength of fluid diffusion within pores to characterize the impact of burial depth on storage. The diffusion coefficient is a key physical quantity commonly used to describe the diffusivity of substances. Accordingly, the diffusion coefficients of CO2 within kaolinite pores at different burial depths were calculated using Formula 6 (Aljama and Wilcox, 2017): D = lim t → ∞ r i t − r i 0 2 6 t (6)

Here, r i is the position of atom i at time t , and r i t − r i 0 2 is the mean square displacement (MSD) of fluid molecules at time t , with d d t is the time derivative of 1/6 of the MSD, yielding the diffusion coefficient D .

From Figure 5, it is observed that with increasing burial depth, the mean square displacement of CO₂ gradually increases. Deriving the diffusion coefficients by fitting the curves, Table 4 shows that the diffusion coefficients also increase with burial depth, indicating enhanced diffusion of CO₂ within the pore. This implies more vigorous fluid movement, indirectly suggesting that the increased burial depth leads to frequent collisions between fluid molecules, resulting in desorption phenomena in the previously stable adsorption layers at the surface, consistent with the conclusions reflected in Figure 6. Moreover, Table 4 reveals that the most significant increase in diffusion coefficient (16.7%) occurs when the burial depth increases from 1.2 km to 2.4 km, indicating a substantial improvement in CO₂ diffusion within the pore when the burial depth exceeds 2.4 km.

To investigate the flow behavior of CO₂ fluid within pores, non-equilibrium molecular dynamics (NEMD) was employed to simulate the experimental pressure drop conditions by applying a fixed acceleration to the horizontal direction of the carbon atoms in each CO₂ molecule. This study explores the influence of different burial depths on the distribution and flow behavior of CO₂ fluid within kaolinite pores.

Firstly, the density distribution of CO₂ at different burial depths within kaolinite pores under a constant pressure drop was explored. Figure 6 reveals that, as CO₂ fluid flows near the surface, stable double adsorption peaks are still formed. With increasing burial depth, the first adsorption peak of CO₂ fluid experiences a slight reduction, decreasing from 1.294 g/cm3 to 1.167 g/cm3, while the density of the second adsorption peak and the bulk phase remains largely unchanged. Applying a constant pressure drop intensifies the motion of each CO₂ molecule, increasing the fluid velocity within the channel, enhancing molecular collisions, and promoting contact between adsorbed and bulk phase molecules. Although the fluid movement within the pore is relatively active, it is insufficient to disrupt the van der Waals forces within the adsorption layers. Consequently, only a small number of CO₂ molecules desorb from the adsorption layers. Figure 7 further indicates that the impact of burial depth on the flowing fluid is minimal, as the adsorption layer density near the surface experiences a slight decrease but stabilizes overall. This suggests that the “effective” storage volume of flowing fluid at different burial depths remains essentially constant, contributing to improved accuracy in estimating the total geological storage capacity of CO₂.

To better elucidate the flow behavior of fluid molecules within pores, this section investigates the velocity distribution of carbon dioxide (CO₂) molecules at different burial depths after applying a constant pressure drop. As shown in Figure 9, the overall velocity profile exhibits a “parabolic” shape, in accordance with the fundamental principles of classical Poiseuille flow, consistent with findings by Wang et al. (Wang and Farzam, 2016). When CO2 fluid passes the kaolinite surface, the restraining effect exerted by the surface diminishes with increasing distance from the surface. Consequently, the fluid velocity gradually increases, reaching a maximum in the central region of the pore. The fitted curve of the distribution graph exhibits a “parabolic” characteristic. Additionally, the velocity of the fluid near the kaolinite surface is greater than zero, indicating a significant slip phenomenon. Despite the strong van der Waals forces between CO₂ and the kaolinite surface, the applied driving acceleration overcomes this binding effect, resulting in fluid slip at the surface. Figure 8 also reveals that with increasing burial depth, the overall fluid velocity increases, showing a subtle upward trend in the velocity curve. At a burial depth of 4.8 km, the maximum velocity is only 9.9% higher than that at a depth of 1.2 km. This suggests that while burial depth promotes fluid velocity, the fluid transport speed primarily depends on the initial pressure drop conditions, consistent with the density distribution patterns observed in Figure 6.

To better elucidate the flow behavior of fluid molecules within pores, this section investigates the velocity distribution of carbon dioxide (CO₂) molecules at different burial depths after applying a constant pressure drop. As shown in Figure 8, the overall velocity profile exhibits a “parabolic” shape, in accordance with the fundamental principles of classical Poiseuille flow, consistent with findings by Wang et al. (Wang and Farzam, 2016). When CO₂ fluid passes the kaolinite surface, the restraining effect exerted by the surface diminishes with increasing distance from the surface. Consequently, the fluid velocity gradually increases, reaching a maximum in the central region of the pore. The fitted curve of the distribution graph exhibits a “parabolic” characteristic. Additionally, the velocity of the fluid near the kaolinite surface is greater than zero, indicating a significant slip phenomenon. Despite the strong van der Waals forces between CO₂ and the kaolinite surface, the applied driving acceleration overcomes this binding effect, resulting in fluid slip at the surface. Figure 8 also reveals that with increasing burial depth, the overall fluid velocity increases, showing a subtle upward trend in the velocity curve. At a burial depth of 4.8 km, the maximum velocity is only 9.9% higher than that at a depth of 1.2 km. This suggests that while burial depth promotes fluid velocity, the fluid transport speed primarily depends on the initial pressure drop conditions, consistent with the density distribution patterns observed in Figure 7.

To accurately describe the slip phenomenon in fluid flow with varying burial depth, this study utilizes the slip length (Ls) to characterize the strength of slip at the wall. Its definition (Liao et al., 2019) is given by Formula 7: L s = v s l i p ∂ v s l i p / ∂ z (7)

Where, v s l i p is the slip velocity of the fluid at the wall, and ∂ v s l i p / ∂ z is the velocity gradient in the z-direction at the wall.

As shown in Figure 10, after applying a constant pressure drop, slip phenomena occur at the pore walls for all fluid locations, indicating that despite the strong van der Waals forces between CO₂ and the kaolinite surface, the applied driving acceleration is sufficient to overcome this binding effect, resulting in slip at the wall. With increasing burial depth, the overall fluid velocity in space increases, and the CO₂ velocity near the wall also expands, leading to a noticeable increase in slip length. Especially at burial depths of 3.6 km and 4.8 km, the slip length significantly increases, aligning with the conclusion of higher diffusion coefficients at greater burial depths as presented in Table 5.