Modeling Adsorption in Silica Pores via Minkowski Functionals and Molecular Electrostatic Moments

Abstract

Capillary condensation phenomena are important in various technological and environmental processes. Using molecular simulations, we study the confined phase behavior of fluids relevant to carbon sequestration and shale gas production. As a first step toward translating information from the molecular to the pore scale, we express the thermodynamic potential and excess adsorption of methane, nitrogen, carbon dioxide, and water in terms of the pore’s geometric properties via Minkowski functionals. This mathematical reconstruction agrees very well with molecular simulations data. Our results show that the fluid molecular electrostatic moments are positively correlated with the number of adsorption layers in the pore. Moreover, stronger electrostatic moments lead to adsorption at lower pressures. These findings can be applied to improve pore-scale thermodynamic and transport models.

1. Introduction

Phase behavior of fluids under confinement is a key process for many natural phenomena and engineering applications [1,2,3,4,5,6,7]. In porous media, this behavior is driven by capillary processes that control the transport and thermodynamic properties of fluids. One example is capillary condensation, the transition of a gas under tight confinement into a condensed phase at a pressure lower than the bulk vapor pressure [8]. This phenomenon is important for estimating hydrocarbon reserves in unconventional reservoirs [1,9] and understanding geologic carbon storage [1,2]. However, translating the impacts of capillary condensation across length- and timescales remains an open area of research [10]. The continuum-scale models become invalid in micropores, and molecular simulations become too expensive to be applied to entire pore networks. In this paper, we explore how integral geometry can effectively be applied as a scale translation approach. We show that a mathematical framework based on Minkowski functionals, using information about the pore morphology, can accurately describe the adsorption behavior of real fluids under tight confinement.

The mechanism behind capillary condensation has been the subject of numerous investigations [1,11,12,13]. This mechanism is strongly influenced by both fluid–fluid interactions and fluid–solid interactions. First, a layer of adsorbed gas forms on the pore walls as a second-order phase transition. In order for this layer of adsorbed gas to transition into a pore-filling condensed phase, the fluid needs to form a liquid bridge in the pore. This formation requires overcoming a nucleation barrier before capillary condensation proceeds as a first-order phase transition. Desgranges and Delhommelle [13] explored the nucleation of liquid bridges and vapor bubbles in confined carbon dioxide (CO₂) and discovered a negative relationship between the chemical potential of coexistence and the fluid–wall interaction. Due to the high density of the condensate, once capillary condensation occurs, the pore can store significantly more fluid molecules. Chen et al. [9] suggested that unconventional gas reserves could be up to six times larger than previously believed because of capillary condensation. In an attempt to better represent phase behavior under confinement, Yang et al. [10] proposed a modified Kelvin equation that reproduces experimental data for the suppressed capillary condensation pressure of various fluids, including CO₂ and hydrocarbon mixtures. Their decomposition showed that the most significant contribution in the shift comes from the fluid–wall interaction. Singh and Mehta [14] studied how two models for the shale solid matrix—graphite and mica—influence capillary condensation in n-alkanes. They showed that mica confinement shrinks the phase envelope of alkanes more than graphite confinement.

The shale solid matrix is heterogeneous and its composition varies widely from one formation to another. Clay minerals, kerogen, and quartz comprise the majority of tight shale. Quartz has been reported to constitute between 30% and 70% of the solid matter of shale [15,16]. Therefore, it is often employed to model shale rock in both experimental and computational studies [1,11,12,17,18,19,20,21]. Capillary condensation—and more broadly adsorption behavior—in quartz pores has been explored for multiple fluids, including water (H₂O) [22,23,24], CO₂ [12,13,20,21], methane (CH₄) [20,21], and nitrogen (N₂) [17,18,25], as common gases found in unconventional gas formations. Burgess et al. [12] reported on adsorption/desorption hysteresis loops for CO₂ in Vycor and connected the hysteresis critical temperature to capillary condensation theory. They also showed that confinement shrinks the phase envelope. Morishige et al. [17,18] employed the volumetric approach to study gas adsorption, including N₂ and CO₂, in silica-based materials. They showed that the hysteresis loop of capillary condensation disappears at higher temperatures, but the point where this occurs is different from the pore critical temperature of the fluid. Raju et al. [26] computationally studied the phase transition of water in graphene nanocapillaries and carbon nanotubes. They related multilayer formation to hydrogen bonding patterns and identified multiple phase transition regimes leading to different crystal structures of confined ice. Kim et al. [27] employed atomic force microscopy to probe the validity of Kelvin’s equation for water between a quartz microscope tip and a mica surface, finding that this thermodynamically derived expression holds for nanopores if the curvature dependence of the surface tension is appropriately accounted for. Water in silica nanopores was also examined by Bonnaud et al. [24]. These researchers, via grand canonical Monte Carlo (GCMC) simulations, observed the adsorption/desorption hysteresis loop and its collapse with smaller pore size. They also employed molecular dynamics (MD) simulations to show that the self-diffusion coefficient of water decreases under confinement, because water molecules adopt an H-down orientation due to the hydrophilic tendency of silica. This layering of the adsorbate in silica pores has been reported for other fluids too. In particular, Jin and Firoozabadi [28] investigated the composition of adsorption layers in clay pores. Studying mixtures of CH₄, CO₂, and H₂O, they reported water forms the layer closest to the pore wall, CO₂ can form multiple adsorption layers (of diminishing adsorption stability), and CH₄ fills in the middle of the pore. Similar behavior was observed by Sun et al. [20] in quartz pores, where they studied the effects of pressure and surface functionaliztion on competitive adsorption of CH₄ and CO₂. Aljamaan et al. [15] experimentally observered multilayer sorption of CO₂ in shale, and suggested that the number of layers depends on the quadrupole moment of the molecule.

The disordered porous nature of the rock substantially impacts the thermophysical properties of confined fluids, introducing large deviation from the properties of bulk fluids. A number of mathematical models have been developed to describe phase behavior under confinement [25,29,30]. Gommes and Roberts [31] proposed a multiphase extension to the standard Gaussian random field model of disordered materials and used this new model to explore capillary condensation. Their findings show that metastable states of adsorption are only found in the most ordered systems considered in that study. Importantly, disorder acts as an environmental noise source. Spagnolo et al. [32] demonstrated that environmental noise affects metastable states in condensed matter systems. These findings are significant for understanding confined fluid behavior in real geologic formations at larger scales; however, the adsorption physics is still unknown at the molecular level. In our work, cylindrical pores are used as a stepping stone in understanding the effects of fluid–solid interactions. The extensions and limitations of our model are also discussed.

Capillary condensation can also be described via Minkowski functionals [29], a mathematical concept for characterizing the morphology and topology of spatial patterns [33], including porous materials. Early studies on this issue have focused on two-dimensional (2D) problems. For a 2D system, there are three Minkowski functionals, representing area, circumference, and signed curvature. Mecke and Arns [34] applied Minkowski functionals to study capillary condensation in a system consisting of overlapping spheres and ellipsoids, successfully reproducing the shift of the critical point under confinement. Boelens and Tchelepi [25] showed that the grand potential (also known as the Landau free energy [35]) and the excess adsorption can be expressed via empirically-augmented formulas based on Minkowski functionals. They adapted Hadwiger’s theorem and classical density functional theory (DFT) calculations with a Lennard-Jones potential. Employing the Minkowski functionals mathematical formalism, they successfully reproduced the phase behavior of nitrogen in 2D smooth Vycor pores.

We extend the work of Boelens and Tchelepi [25] and examine the validity of their approach for more complex interatomic potentials, rough surfaces, and four different compounds with different molecular size and polarity. We show that the Minkowski functionals can accurately describe the phase behavior of real fluids in SiO₂ nanopores despite the competing effects of the various molecular interactions. Our results demonstrate that Minkowski functionals can be used to translate information from the molecular level to pore scale. In exploring the impact of molecular size and polarity on adsorption, we find that the adsorption behavior of a fluid is strongly influenced by the dipole and quadrupole moments of the fluid molecules. The dipole moment is a measure of the polarity of a chemical bond. The quadrupole moment describes the effective shape of the molecular charge distribution. A perfectly spherically symmetric charge distribution has zero quadrupole moment. In this work, we show that stronger electrostatic moments lead to early onset of adsorption and to the formation of multiple adsorption layers across the pore.

The remainder of this paper is organized as follows. Section 2 discusses how Minkowski functionals can be applied to three-dimensional (3D) systems and outlines our molecular simulations. Section 3 presents the findings of our study and discusses their scientific significance. Section 4 summarizes the main results and discussed their implications for geologic carbon storage and gas production from unconventional formations.

2. Methods

2.1. Mathematical Framework

Following the approach of Boelens and Tchelepi [25], we aim to express the grand potential in terms of Minkowski functionals. The grand potential $\Omega =\Omega [X;\mu ,T]$ is an extensive quantity that depends on the geometric properties of the system X, the chemical potential $\mu$, and the temperature T. Thermodynamically, we can completely describe a system if we know its appropriate thermodynamic potential, that is, the grand potential in the case of the grand canonical ensemble. By separating variables and assuming linearity, the grand potential can be written as:

$$\Omega [X;\mu ,T]=\omega (\mu ,T)V\left[X\right],$$

(1)

where $\omega (\mu ,T)$ is an intensive thermodynamic quantity and $V\left[X\right]$ is the volume of the system. For the grand potential, $\omega (\mu ,T)$ is the negative of the pressure $p(\mu ,T)$ of the system. Equation (1) is valid only for bulk systems. For systems under confinement, the dependence of $\Omega [X;\mu ,T]$ on the shape of the system X is an integral over the phase space of the system. For real pores, this integral can only be approximated [36]. To overcome this limitation, König et al. [36] suggested three physical restrictions on this dependence (shown in Figure 1):

• Motion invariance: The thermodynamic potential is independent of the system’s location and orientation in space. For any translation or rotation g:

  $$\Omega \left[X\right]=\Omega \left[gX\right].$$

  (2)
• Continuity: If a sequence of convex sets $X_n$ converges to X for $n\to \infty$, then $\Omega \left[X_n\right]\to \Omega \left[X\right]$. This property states that an approximation of the convex domain also yields an approximation of the grand potential.
• Additivity: The functional union of two domains $S_1$ and $S_2$ is the sum of the functional of the individual domains subtracted by their intersection:

  $$\Omega [S_1\cup S_2]=\Omega \left[S_1\right]+\Omega \left[S_2\right]-\Omega [S_1\cap S_2].$$

  (3)

Enforcing these properties, one can invoke Hadwiger’s theorem [34,36], stating that any motion-invariant, continuous, and additive functional can be written as a linear combination of the corresponding Minkowski functionals. Specifically, for a 3D system the Minkowski functionals are defined as [29]:

$$\begin{array}{ccc}\hfill & M_0\left[X\right]=V\left[X\right]\hfill & \end{array}$$

(4)

$$\begin{array}{ccc}\hfill & M_1\left[X\right]=\frac{1}{8}A\left[X\right]\hfill & \end{array}$$

(5)

$$\begin{array}{ccc}\hfill & M_2\left[X\right]=\frac{1}{2{\pi }^2}H\left[X\right]\hfill & \end{array}$$

(6)

$$\begin{array}{ccc}\hfill & M_3\left[X\right]=\frac{3}{4\pi }\chi \left[X\right],\hfill & \hfill \end{array}$$

(7)

where V is the volume of a pore space X, A is the surface area of the pore, H is the integral mean curvature, and $\chi$ is the Euler characteristic of the pore. Using Hadwiger’s theorem, Mecke et al. [34,36] give the grand potential as:

$$\Omega =-p(\mu ,T)M_0\left[X\right]+\sigma (\mu ,T)M_1\left[X\right]+\kappa (\mu ,T)M_2\left[X\right]+\overline{\kappa }(\mu ,T)M_3\left[X\right],$$

(8)

where $\sigma (\mu ,T)$ is the surface tension and $\kappa (\mu ,T)$ and $\overline{\kappa }(\mu ,T)$ are bending rigidities. These coefficients are properties of the fluid and fluid–wall interactions, but are independent of the geometry of the wall.

We note that for an infinite cylinder, the integral mean curvature $H=1/\left(2R\right)$, where R is the pore radius, and the Euler characteristic is equal to zero. Therefore, we use the following expressions for the Minkowski functionals of the cylindrical pores in our study: $M_0\left[X\right]=\pi R^{2}L$, $M_1\left[X\right]=\pi RL/4$, $M_2\left[X\right]=1/\left(4{\pi }^{2}R\right)$, and $M_3\left[X\right]=0$. Thus, the last term in Equation (8) is zero for the geometries investigated in our work. In these expressions, L is the length of the simulation box along the pore central axis. Based on numerical results, Boelens and Tchelepi [25] augmented this expression for Lennard-Jones fluids by adding two additional terms. One of these terms is proportional to $\sqrt{M_0\left[X\right]}$ and the other term is proportional to $M_1^{3/4}\left[X\right]$. Due to the non-linearity, these two terms are not additive when studying pore networks [25,29]. In contrast to the 2D DFT simulations of Boelens and Tchelepi [25], our GCMC approach models the pore as a cylinder in 3D. Therefore, we translate their empirical terms to 3D and rewrite the augmented version of Equation (8) as:

$$\Omega =-p(\mu ,T)M_0\left[X\right]+\sigma (\mu ,T)M_1\left[X\right]+\kappa (\mu ,T)M_2\left[X\right]+p^{\prime }(\mu ,T)M_0^{1/2}\left[X\right]+{\sigma }^{\prime }(\mu ,T)M_1^{3/4}\left[X\right].$$

(9)

The coefficients of the last two terms are named pseudo-pressure and pseudo-surface tension [25]. All coefficients in this expression are functions of the chemical potential and temperature, but not the geometry. Therefore, for each fluid, there is a single set of coefficients independent of pore size. These coefficients are obtained, following the approach of Boelens and Tchelepi [25], by performing a least-squares fit for all our numerical data to Equation (9).

All relevant thermodynamic properties in the grand canonical ensemble can be computed from the grand potential [37]. In laboratory experiments, for example, a more commonly measured quantity is the excess adsorption ${\Gamma }_{\mathrm{ex}}$. To calculate this quantity, we first use the grand potential to find the interfacial tension:

$$\gamma =\frac{\Omega +p_{B}V}{A},$$

(10)

as a measure of the change in the grand potential introduced by the wall per unit area. Here, $p_B$ is the pressure in the bulk phase. In terms of Minkowski functionals, we can express $\gamma$ as [25]:

$$\gamma =(p_B(\mu ,T)-p(\mu ,T))\frac{M_0\left[X\right]}{M_1\left[X\right]}+\sigma (\mu ,T)+\kappa (\mu ,T)\frac{M_2\left[X\right]}{M_1\left[X\right]}+p^{\prime }(\mu ,T)\frac{M_0^{1/2}\left[X\right]}{M_1\left[X\right]}+{\sigma }^{\prime }(\mu ,T)M_1^{-1/4}\left[X\right].$$

(11)

Finally, using Gibbs theorem, we obtain the excess adsorption [25]:

$${\Gamma }_{\mathrm{ex}}=-{\left(\frac{\partial \gamma }{\partial \mu }\right)}_{T,V}=\frac{\partial }{\partial \mu }(p-p_B)\frac{M_0\left[X\right]}{M_1\left[X\right]}-\frac{\partial \sigma }{\partial \mu }-\frac{\partial \kappa }{\partial \mu }\frac{M_2\left[X\right]}{M_1\left[X\right]}-\frac{\partial p^{\prime }}{\partial \mu }\frac{M_0^{1/2}\left[X\right]}{M_1\left[X\right]}-\frac{\partial {\sigma }^{\prime }}{\partial \mu }{M}_1^{-1/4}\left[X\right],$$

(12)

where the coefficients in front of the Minkowski functionals are the derivatives with respect to chemical potential of the coefficients in Equation (9). In this paper, we explore how effectively this analytical model captures observations from molecular simulation.

2.2. Numerical Simulations

The numerical results presented in this work are obtained with the Large-Scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) software package [38,39]. All simulations are performed in the grand-canonical ($\mu VT$) ensemble and are based on Monte Carlo exchanges with an ideal gas reservoir of the fluid at the specified temperature and chemical potential. The temperature is controlled using a Nosé-Hoover thermostat with a damping constant of 50 fs. The temperature of the simulations is $T/T_c=0.61$, where $T_c$ is the bulk critical temperature of the fluid. This value was chosen based on the convention of studying adsorption with nitrogen at 77 K [40]. In Kelvin, the simulations’ temperature is 116.8, 77, 185.3, and 389.5 K, for simulations involving methane, nitrogen, carbon dioxide, and water, respectively. For each of the four fluids, we vary the chemical potential between $\mu /k_{B}T=-200$ and $-1$. A non-dimensional interval of $\Delta \mu /k_{B}T=1$, dimensionally is 0.23, 0.15, 0.37, 0.77 kcal/mol for methane, nitrogen, carbon dioxide, and water, respectively.

The bulk simulations are done in a simulation box of $20\times 20\times 20$ Å${}^3$. The confined simulations are done in pores with diameters d = {5, 10, 15, 20, 40} Å. Each pore is 31.8 Å long. In both bulk and confined simulations, periodic boundary conditions are used. The number of particles in the pores varies between 0 and 1800 molecules. All simulations are equilibrated for at least 400,000 steps, and statistics are collected over 200,000 steps.

The interactions between particles are represented using the TraPPE-UA force field for methane [41], TraPPE for carbon dioxide and nitrogen [42], SPC/E for water [43], and ClayFF for the SiO₂ pores [44]. TraPPE treats CO₂ and N₂ as rigid bodies, while the bond distances and bond angles of SPC/E H₂O are constrained via the SHAKE algorithm. TraPPE-UA CH₄ is represented via the unified-atom approach. The force field parameters used are summarized in

Table 1. Summary of the force field parameters used in this work. References to the force fields are given in the text.

Molecule	Particle Type	$\mathit{\epsilon }/{\mathit{k}}_{\mathit{B}}$ (K)	$\mathit{\sigma }$ (Å)	q (e)	${\mathit{k}}_{\mathit{bond}}$ (kcal/mol)	${\mathit{k}}_{\mathit{\theta }}$ (kcal/mol)
CH4	CH4	148	3.730	-	-	-
CO2	C	27	2.800	0.7000	5000.0	500.0
	O	79	3.050	−0.3500
H2O	H	0.0	0.000	0.4236	1000.0	100.0
	O	78.2	3.166	−0.8472
N2	N	36	3.310	−0.4820	5000.0	500.0
	Ghost	0.0	0.000	0.9640
SiO2	Si	$9.26\times {10}^{-4}$	3.302	2.1000	-	-
	O	78.2	3.166	−1.0500	-	-

. The non-bonded interactions are given by the Lennard-Jones potential and Coulombic electrostatic interactions. A 12 Å cut-off radius is used for both potentials. The long-range electrostatic interactions are calculated with the Ewald summation method with relative tolerance of ${10}^{-4}$ for forces. To validate that these potentials are appropriate for studying adsorption, in Figure 2, we compare the adsorption isotherm of N₂ at 77 K in a 20-Å-diameter silica pore obtained with our GCMC simulations to an experimental measurement by Morishige et al. [17]. There is a good agreement between the two adsorption isotherms. Moreover, similar combinations of force fields have been used previously for studying the transport and thermodynamics of confined fluids [45,46,47,48], further validating the choice of potentials.

3. Results and Discussion

We study the phase behavior under tight confinement of four fluids (CH₄, N₂, CO₂, and H₂O) with different electrostatic moments, which are summarized in

Table 2. Summary of the molecular electrostatic moments for investigated fluids. The dipole is given in Debye (D) and the quadrupole moment is in Debye-Ångstrom (DÅ).

Molecule	Force Field	Electrostatic Dipole (D)	Electrostatic Quadrupole (DÅ)	Reference
CH4	TraPPE-UA	0.0	0.0	-
N2	TraPPE	0.0	−1.47	[49]
CO2	TraPPE	0.0	−4.1	[50]
H2O	SPC/E	1.86	0.11	[51]

. These fluids commonly occur in shale formations. For our study, the fluids are confined in cylindrical SiO₂ pores with diameters ranging between 5 and 40 Å, representative of tight shale rocks. Analyzing the molecular trajectories of GCMC simulations, we show that the onset of capillary condensation is uniquely defined by pore morphology. As shown in Figure 3, fluid molecules form one or multiple adsorption layers on the pore surface before filling up the entire pore. The number of layers formed depends on the electrostatic moments of the molecule. Finally, we connect the thermodynamics and pore morphology by extending the Minkowski functionals expression proposed by Boelens and Tchelepi [25] to a wider range of fluids and to 3D systems.

3.1. Influence of Molecular Polarity on Adsorption

Figure 3, Figure 4 and Figure 5 show that molecules with large electrostatic moments form more adsorption layers near the pore surface. For each fluid, Figure 4 shows a map of the observed phase behavior in GCMC simulations. In particular, the region where multi-layer adsorption occurs broadens going from methane (no molecular electrostatic moments) to water (polar molecule). Larger molecular electrostatic moments lead to stronger interaction between the silica pore and the fluid, which in turn leads to a more stable physisorption of fluid molecules on the pore walls. Generally, in each map we observe four regions: an empty pore, adsorption in a single layer near the pore wall, forming multiple adsorption layers near the pore wall, and an entirely filled pore. On one hand, the region of entirely filled pore does not evolve significantly with the polarity of different species. On the other hand, the regions of single-layer adsorption and multi-layer adsorption grow considerably for increasing molecular electrostatic moments. These two regions are very limited for methane, where we typically see immediate filling of the entire pore (second row of Figure 3); yet, they span a broad range of chemical potential values for water, where we see gradual formation of adsorption layers before the pore is entirely filled.

Moreover, as shown in Figure 5, the density in the interlayer regions increases for fluids with larger molecular electrostatic moments. This increase implies that there is close packing of the molecules even between the defined adsorption layers. The first row of snapshots in Figure 3 shows that CO₂ fills the pore with closely packed layers. This is in contrast to CH₄, whose adsorption layers can be easily delineated in the second row of Figure 3.

The effect of the electrostatic moments on the close packing of the adsorption layer is countered by the size of the pore and the molecule. This is illustrated in Figure 5, showing wider valleys in the transverse density profiles for fluids with larger molecular diameters. The transverse density profiles of CO₂ and N₂, whose molecular diameters differ by about 10%, show similar structure with wider valleys for CO₂. Smaller pore sizes also impose a limitation on the adsorption behavior. The stronger, overlapping fluid–wall interactions lead to lower average fluid density in small pores. This dependence of the fluid density on pore size beyond the onset of capillary condensation suggests decreased accessible free volume because of repulsive wall effects and geometric constraints.

3.2. Grand Potential

Following Peterson and Gubbins [52], we determine the grant potential from the GCMC simulations by integrating the Gibbs adsorption isotherm:

$${\left(\frac{\partial \Omega }{\partial \mu }\right)}_{T,V}=-N,$$

(13)

where N is the average number of molecules in the pore. The initial condition used in solving this ordinary differential equation (ODE) is the grand potential of an ideal gas [37]. We solve the ODE using an explicit Runge-Kutta integrator in MATLAB [53]. The results are shown as solid lines in Figure 6. We start studying the adsorption at very low chemical potentials when all pores are empty. As we increase the chemical potential, the slope of the grand potential becomes steeper, because molecules start adsorbing in the pore. In the range between $\mu /k_{B}T=-25$ to $-15$, depending on the fluid, the slope changes sharply, indicating a jump in the density inside the pore and onset of the capillary condensation as a first-order phase transition.

Using Equation (9), we determine the grand potential in terms of Minkowski functionals. In Figure 6, we compare the numerical results from GCMC simulations (solid lines) with this geometric reconstruction (symbols). We see excellent agreement for the four fluids in all pore sizes. The geometry-based reconstruction is able to capture phase behavior of both polar and non-polar fluids in very tight confinement. As the electrostatic moments of the fluid molecules change (Table 2), their phase behavior in the pore becomes more complex (Figure 4). Methane exhibits a single discontinuous phase transition, accompanied by a sharp change in the slope of the grand potential. Nitrogen shows three distinct regimes in its behavior: an empty pore (zero slope for $\mu /k_{B}T<-90$), a layer of adsorbed gas on the pore surface (slowly varying slope for $-90<\mu /k_{B}T<-18$), and condensed fluid filling the pore (constant steeper slope for $-18<\mu /k_{B}T$). These three regimes appear as three different slopes in the grand potential. Carbon dioxide exhibits multiple phase regimes as layers of gas form near the pore surface until the pore is filled (see Figure 3). At that point, the grand potential obtains a constant slope. Finally, water shows continuously varying slope with a small discontinuous jump in fluid density (change in slope of the grand potential) at high pressures (i.e., high chemical potential). This jump is the onset of capillary condensation. The Coulombic interactions between the fluids and the pore walls define the molecular orientation within the layers and allow for closer packing of the fluid. Water as a polar molecule with a dipole moment forms more layers than methane, nitrogen, and carbon dioxide (see Figure 5).

This close packing in the adsorption layers is counteracted by the small pore sizes. In the range of length scales explored in this work, the pore diameter, which is on the order of few molecular diameters, physically constrains the number of layers that can be formed. Therefore, in smaller pores, we see capillary condensation occurring at lower pressure (i.e., lower chemical potential), accompanied by an earlier discontinuous change in the slope of the grand potential.

3.3. Excess Adsorption

The excess adsorption is a physical quantity that defines the additional amount of fluid adsorbed in a pore of given volume as compared to a bulk fluid at the same thermodynamic condition [54]. We normalize this quantity with the surface area of the pore, such that the impact of tight confinement can be discussed. To calculate the excess adsorption from GCMC simulations, we run both simulations of confined fluid and of bulk fluid and use those density values in the following equation:

$${\Gamma }_{ex}=\frac{1}{A\left(X\right)}\left(\rho (\mu ,T)-{\rho }_B(\mu ,T)\right)V\left(X\right).$$

(14)

Figure 7 compares the excess adsorption obtained from numerical simulations (solid lines) with its geometric reconstruction using Minkowski functionals (symbols). The two compare very well for all four fluids. The results show that the pore size directly affects the extent to which fluid can be adsorbed, mostly by preventing fluid molecules from entering the smallest pores because of their physical size. However, in pores of 10-Å diameter, the strong fluid–wall interactions push most molecules to the pore center and limit the formation of well-defined adsorption layers (Figure 5). These interactions result in less excess adsorption in pores whose geometrical size does not mechanically act as a molecular sieve. Thus, in the range of conditions where only adsorption layers exist in the pore, the storage capacity does not substantially change with pore size. But, we emphasize that because of the larger volume-to-area ratio, larger pores can store more fluid after the onset of capillary condensation (Figure 7), supporting the enhanced hydrocarbon estimates by Chen et al. [9].

Fluids with large electrostatic moments, such as water, start to adsorb at significantly lower pressure (i.e., lower chemical potential) and form more adsorption layers before the onset of capillary condensation than fluids with no or very low electrostatic moments, such as methane. To further demonstrate that the electrostatic interactions lead to this difference in adsorption behavior, we compare the adsorption behavior of two nitrogen models. These additional GCMC simulations model nitrogen without atomic charges. This model was used in the DFT simulations of Boelens and Tchelepi [25] and follows the unified-atom (UA) approach. The particles interact solely via the Lennard-Jones potential with $\epsilon /k_B=94.45$ K and $\sigma =3.575$ Å. Figure 8A compares the bulk fluid density of TraPPE N₂ (used in our study) and the UA N₂ model. Both models accurately capture the phase transition of bulk nitrogen at 77 K. Next, Figure 8B shows the fluid density inside a 10-Å-diameter SiO₂ pore for both models. The significant difference between the density curves in the pore suggests that electrostatic forces have a strong influence on adsorption for non-polar fluids. On one hand, TraPPE nitrogen, where individual atoms have charge and the molecule has a quadrupole moment, starts adsorbing at low pressures and forms adsorption layers as the pressure increases. On the other hand, the UA nitrogen model sees minimal adsorbed fluid amount before a sharp phase transition at approximately $\mu /k_{B}T=-22$, when the pore fills up. Therefore, we conclude that any proposed framework for describing phase behavior in nanoporous materials must accurately capture the impact of the molecular electrostatic interactions. Similar observations were reported by Hackett and Hammond [55], who examined the adsorption of different nitrogen models in zeolites.

3.4. Toward Real Systems

This work has focused on cylindrical SiO₂ pores as a model for quartz grains in the shale solid matrix. Similar models have been informative for reproducing observations from real shale rocks with similar pore size and geochemical composition [56]. Yet, the nature of the shale pore network is complex and disordered [11]. The variation of pore morphology and topology in the shale rock influences the conditions under which capillary condensation would occur in the system as a whole. A number of models have been developed to explore this relationship in a statistical manner for heterogeneous, disordered systems at the meso- and larger length scales [30,31,57]. Despite the significance of these theoretical findings, Gommes and Roberts [31] call for further validation of their model against experimental small-angle scattering data. The disorder of the porous material acts as an environmental noise source, which has been reported to have a positive role in condensed matter systems. Spagnolo et al. [32] found that noise can extend the lifetime of metastable states in fluctuating systems. Future studies are needed to explore how environmental noise sources affect metastable states of sorption phenomena under tight confinement. Coasne et al. [11] discussed the phase behavior of fluids confined in porous silica and considered how morphological and topological disorder affects phase phenomena. These effects include delayed phase transitions and asymmetric sorption hysteresis loops. However, the authors acknowledge the significance of characterizing confined fluid behavior in regular and ordered systems as a crucial step in developing more complex, realistic models. Therefore, describing capillary condensation in cylindrical pores and other simplified geometries is an important building block for understanding confined-fluid phenomena in shale and other geologic formations [56,58,59]. The Minkowski functionals framework presented in this work can be used in the future to calculate excess adsorption and the grand potential in disordered pore networks. The Minkowski functionals are well defined (Equations (4)–(7)) for all 3D geometric shapes. Thus, the only challenge with applying this framework to a disordered pore network is theoretically resolving the additivity of the non-linear terms in Equations (9) and (12).

Furthermore, the molecular approach to studying fluid phase behavior is very sensitive to environmental fluctuations in the system. The molecular velocities in equilibrium systems always follow a distribution at finite temperatures. Experimentally, these fluctuations are enhanced by limitations in the control of macroscopic thermodynamic variables. In GCMC simulations, the computational thermostat also marks fluctuations due to the random forces [60]. Moreover, the number of molecules in the simulation box fluctuates around a mean value, introducing density fluctuations. Similar fluctuations in the number of molecules inside the pore occur in adsorption experiments. Tanaka et al. [61] argued that fluctuations in experiments are larger than those in GCMC simulations, thus affecting the observed phase behavior. These fluctuations, in density or energy of the system, become large near the critical point. Furthermore, they can lead to a phase transition by pushing the system over the energy barrier for nucleation of a condensed liquid bridge or a vapor bubble inside the pore. To ensure that our simulations accurately represent the phase behavior in silica pores, we varied the chemical potential in smaller increments in regions with large change in number of molecules. Moreover, we allowed for sufficient time for the system to equilibrate inside the pore and collected statistics over 200,000 steps, as discussed in Section 2.

When studying confined fluid behavior, another important consideration is the difference in thermophysical properties between confined and bulk fluids. This difference is a result of the complex fluid–solid interactions that occur inside the pores, especially under tight confinement, in addition to the fluid–fluid interactions that are found in bulk fluids. One example of how confinement influences fluids’ thermophysical properties is the shift in the phase envelope. This shift has been well-documented for various fluids in the literature [12,14]. While it is out of scope for the present study to map the shift of the phase envelope across a range of temperatures, we observe the gas-liquid transition occurring at different chemical potentials (different pressures) along the studied isotherms as a function of the pore diameter. Most notably, the smallest pores undergo first-order liquid-gas transition first at a low chemical potential. As the pore size is increased, the chemical potential at which a gas-liquid phase transition occurs is monotonously growing (Figure 7).

The density of the confined liquid phase fluid is also different from the liquid phase fluid at the same chemical potential and temperature in the bulk. Chiang et al. [58] employed small-angle neutron scattering measurements to show that the density of liquid deuterated methane (CD₄) in mesoporous silica is lower than that of liquid bulk CD₄. A similar conclusion was reached theoretically by Tan and Piri [62], who studied CO₂ and n-pentane under confinement via the perturbed-chain statistical associating fluid theory (PC-SAFT) coupled with the Young-Laplace equation. In our simulations (Figure 7), we report a lower density of confined liquids compared to their bulk equivalents for CH₄, N₂, and CO₂. Water, however, exhibits the opposite behavior where liquid water under confinement is denser than in the bulk at the same temperature and chemical potential. This difference is probably a result of the surface roughness in our model combined with the molecular size. Larger molecules cannot fit within the corrugated pore surface. Moreover, as discussed in the previous section, H₂O exhibits extremely close-packing, while CH₄, CO₂, and N₂ have more pronounced valleys and peaks in their transverse density profiles (Figure 5). A further study is needed to explain why confined liquids are less dense than their bulk equivalents.

4. Conclusions

The adsorption of fluids under tight confinement is crucial for many technological and natural processes, including carbon sequestration and shale gas production. In this work, we examine how four fluids fill SiO₂ nanopores. These fluids—methane, nitrogen, carbon dioxide and water—were chosen because they are commonly found in subsurface reservoirs and aquifers. Moreover, we offer a mathematical framework to describe the phase behavior of confined fluids. Expanding on the work of Boelens and Tchelepi [25], we show that a mathematical framework in terms of Minkowski functionals (volume, surface area, and integral mean curvature) can reproduce capillary condensation with good accuracy. Using geometric information for the porous rock to describe the phase behavior of different fluids paves the way for better scale translation of these processes.

Our work shows that Minkowski functionals can leverage geometric information about the pore to describe thermodynamic quantities such as the grand potential and the excess adsorption. We apply this method to both polar and non-polar fluids, and show that the molecular electrostatic moments are crucial to the adsorption layer formation. Stronger moments lead to early onset of adsorption in the pores and formation of multiple layers before the pore is entirely filled. This observation is significant because it reveals why fluid–solid interactions are important in adsorption phenomena and reasserts the importance of accurately capturing the electrostatic interactions when modeling adsorption phenomena at the molecular level. This paper presents the first application of the Minkowski functionals framework to describe the phase behavior of real fluids in 3D.

The results and the mathematical formulation presented in this work can be extended to more complex geometries. Further research is needed to investigate the effects of pore network disorder [11] and environmental noise [32] on adsorption layer formation under tight confinement. Finally, we reiterate that the non-linear terms of the Minkowski-functionals reconstructions (e.g., last two terms of Equations (9) and (12)) are not additive when a whole pore network is studied [29]. Further theoretical work is needed to address this issue.