Propane-Water Mixtures Confined within Cylindrical Silica Nanopores: Structural and Dynamical Properties Probed by Molecular Dynamics.

Despite the multiple length and time scales over which fluid-mineral interactions occur, interfacial phenomena control the exchange of matter and impact the nature of multiphase flow, as well as the reactivity of C-O-H fluids in geologic systems. In general, the properties of confined fluids, and their influence on porous geologic phenomena are much less well understood compared to those of bulk fluids. We used equilibrium molecular dynamics simulations to study fluid systems composed of propane and water, at different compositions, confined within cylindrical pores of diameter ∼16 Å carved out of amorphous silica. The simulations are conducted within a single cylindrical pore. In the simulated system all the dangling silicon and oxygen atoms were saturated with hydroxyl groups and hydrogen atoms, respectively, yielding a total surface density of 3.8 -OH/nm2. Simulations were performed at 300 K, at different bulk propane pressures, and varying the composition of the system. The structure of the confined fluids was quantified in terms of the molecular distribution of the various molecules within the pore as well as their orientation. This allowed us to quantify the hydrogen bond network and to observe the segregation of propane near the pore center. Transport properties were quantified in terms of the mean square displacement in the direction parallel to the pore axis, which allows us to extract self-diffusion coefficients. The diffusivity of propane in the cylindrical pore was found to depend on pressure, as well as on the amount of water present. It was found that the propane self-diffusion coefficient decreases with increasing water loading because of the formation of water bridges across the silica pores, at sufficiently high water content, which hinder propane transport. The rotational diffusion, the lifespan of hydrogen bonds, and the residence time of water molecules at contact with the silica substrate were quantified from the simulated trajectories using the appropriate autocorrelation functions. The simulations contribute to a better understanding of the molecular phenomena relevant to the behavior of fluids in the subsurface.

Introduction

The physical properties of fluids confined in various environments are of considerable interest, given that they could provide information on glass transition,[1,2] migration[3,4] and adsorption,[5−7] phenomena that are relevant to geology, biology, and engineering. Keith Gubbins and his group pioneered the use of molecular simulations for better understanding the behavior of confined fluids.[8] He addressed many important topics, such as phase separations in confinement,[9] freezing and melting in confinement,[10] and also the development of realistic models for the porous adsorbents.[11] This Article considers fluids in silica nanopores. Several authors have studied structure and dynamics of fluids confined in nanoporous silica materials with different geometries (slit pores[12−16] or cylindrical pores[17−21]). These studies provide insights for applications ranging from nanofiltration,[22] separation,[23] and catalysis,[24] among others.

Micro- and mesoporous silica materials have been widely used for systematic studies of fluids in confined environments, in part because silica is one of the most abundant subsurface materials. When pure silica comes in contact with water, its outer layer is expected to be hydroxylated. The hydrophilic character of the resultant pore surface significantly affects the confined fluids, as it has been documented in the case of interfacial water both in terms of structure and dynamics.[19,20,25−27] The nature of the interactions between water, guest molecules, and the solid substrate must be better understood to elucidate molecular phenomena that occur in narrow pores at subsurface conditions. Toward this goal, our group previously reported molecular dynamics (MD) simulations that documented the effect of density and composition of CO2-hydrocarbon mixtures on adsorption and mobility of the confined mixtures.[4,16] Those simulation results, which seem to be in agreement with experimental observations,[28−30] suggest that the preferential CO2 adsorption to the pore walls weakens the adsorption of hydrocarbons, and enhances the self-diffusion of hydrocarbons via lowering the activation energy for diffusion. The present work stems from recent quasi-elastic neutron scattering experiments conducted for systems containing propane confined in MCM-41-S materials.[31] These experimental results showed that propane mobility decreases as the D2O content increases. To probe a system similar to the experimental one, in this manuscript we consider propane–water mixtures confined inside a 16 Å diameter cylindrical silica pore. The geometry of the pore and its reduced size, the amorphous nature of the pore surface, and the preferential interactions between water and the surface −OH groups as opposed to the weaker interactions between propane and the silica substrate are expected to yield significant differences compared to the results reported previously for CO2-hydrocarbon mixtures in slit-shaped pores carved out of crystalline cristobalite. The simulations were conducted at moderate temperature–pressure conditions (T = 300 K, Pbulk = ∼ 0.6–3 MPa). Our attention focused on the effect of the bulk pressure of C3H8 loading as well as on the effect of adding water on the dynamics of confined C3H8. The simulation results provide extensive insights into the structural and dynamic properties of all components considered.

The remainder of this manuscript is organized as follows: We start by reviewing both simulation models and algorithms, we then present our main simulation results, we compare the results with those from other mixed fluid-porous matrix interactions, and we finally conclude via summarizing our main results.

Simulation Methods and Algorithms

Preparation of the Amorphous Silica Pore

To create a model of bulk amorphous silica, we started from the β-cristobalite structure in a system composed of 12288 atoms within a 57.28 Å × 57.28 Å × 57.28 Å simulation box with periodic boundary conditions. The amorphous system was prepared following the annealing cycle proposed by Leroch et al.[32] In this procedure, the crystalline sample is melted at 7000K, then equilibrated in the liquid phase, and finally quenched to room temperature at a rate of 4K/ps. The annealing simulations were conducted first in the NPT ensemble at a pressure of 1 bar, which allows for adaptations of the simulation cell volume during phase transition to maintain the desired pressure. At the end of annealing process, the amorphous silica block was equilibrated under NVT conditions. The final cubic dimension of the simulation box was of 56.9 Å on all sides.

Following Leroch et al., for the annealing simulations we implemented the Morse-type potential developed by Demiralp et al.[33] In this force field, a two-body Morse–Stretch term describes nonelectrostatic interactions, and the atomic charges depend on the local atomic configuration. The interaction potential is described asIn eq q and q represent the charges of atoms i and j, respectively. In our approach, for Si atoms qSi = +1.3e and for O atoms qO = −0.65e. In eq R denotes the interatomic distance between atoms i and j. D0, R0 and γ represent bond strength, bond length, and dimensionless force, respectively. It has been shown that the potential of eq closely reproduces the melting temperature of cristobalite, the glass phase transition temperature of silica glass, and the density of silica after the annealing cycle,[34−36] in some cases better than BKS[37] or TTAM[38] potentials. In Table we report the Morse potential parameters used in this study.

Table 1

Morse Potential Parameters Taken from Ref (36)

interaction	R0 (Å)	D0 (kcal/mol)	γ
O–O	3.7910	0.5363	10.4112
Si–Si	3.7598	0.17733	15.3744
Si–O	1.6280	45.9970	8.6342

From the bulk amorphous silica volume prepared as described above, a cylindrical channel of 16 Å in diameter was carved out by removing all atoms located within a distance of 8 Å from the X axis, i.e., all atoms whose y and z coordinates obey the relationThe resulting SiO2 substrate was composed of 11613 atoms. The surface of the cylindrical pore is rough at the atomic level, reflecting the amorphous nature of the substrate. The dangling silicon and oxygen atoms in the interfacial region were saturated with hydroxyl groups and hydrogen atoms, respectively, yielding two kinds of silanol groups: Si–(OH)2 (germinal) and Si–OH (single silanol). The resultant hydroxyl density was of 3.8/nm2, which is in good agreement with experimental data (from 2.6 to 4.6 OH/nm2) measured on flat amorphous silica surfaces.[39] The final number of atoms in the solid substrate was 11 694.

Within the periodic boundary conditions implemented, the cylindrical pore was aligned parallel to the X direction, along which it was effectively infinite. This model substrate was used to sample transport and structural properties for the confined fluid systems.

To prepare a model of propane–water mixtures confined inside the pore at equilibrium with a bulk reservoir, we followed the procedure previously implemented for other substrates,[15,40] and briefly described below.

Once the amorphous pores were prepared, subsequent simulations concerning fluid-pore systems were conducted implementing CLAYFF for the solid substrate, as described below. Two systems were prepared and used in our simulations. In one (System A), the solid substrate was in contact with an external “bulk” fluid reservoir. In the other (System B), the pore was infinitely long because of periodic boundary conditions. In Figure we report a schematic of the two simulation scenarios. The two corresponding simulation box sizes were of 150.0 × 56.9 × 56.9 and 56.9 × 56.9 × 56.9 Å3, respectively. In both cases, periodic boundary conditions were implemented along the X, Y, and Z directions. The X dimension of the simulation box for System A simulations allowed for a “bulk” region of thickness ∼93 Å. The X dimension of the simulation box for System B simulations, 56.9 Å, is large enough to prevent unphysical interactions between replicas of the simulated molecules along the X direction. It is worth repeating that we employed System A to determine how much propane is adsorbed within the pore by varying the thermodynamic conditions of bulk propane. The final configuration obtained from simulations of System A was then used to prepare System B, which was employed to quantify the transport and the structural properties of the fluid systems confined in the pore. One alternative approach would be to use grand canonical Monte Carlo simulations to obtain the initial configurations for the simulations conducted within the System B. In a prior work,[6] we showed that the approach followed within the System A algorithm yields results that are in fair agreement with adsorption isotherms.

Figure 1

Side view of representative simulation snapshots for setups following System A (top) and System B (bottom), as described in the text. Yellow, red, white, and cyan spheres represent silicon, oxygen, hydrogen atoms, and propane molecules, respectively. The silica substrate contains 3847 silicon, 7745 oxygen, and 102 hydrogen atoms. Note that only a portion of the simulation box is shown for clarity.

System A

The initial configuration for simulations conducted within this framework consists of a desired number of water and propane molecules placed at each side of the cylindrical pore, in the “bulk” region, along the X direction. As simulations proceeded, water and propane spontaneously filled the pore and were distributed across both pore and bulk volumes. To fill the pore completely with water we require 2141 water molecules. An increasing amount of propane (from 76 to 350 molecules) is placed in the bulk region. This propane controls the pressure of the system. As the simulations progress, some propane adsorbs in the hydrated pore. Once equilibrium is reached, the propane density in the bulk was calculated from density profiles such as those shown in Figure b. These representative results show that, for each simulated system, propane accumulates near the substrate and penetrates the pore, while a constant propane density is maintained away from the solid substrate. The portion of each data set where the propane density is constant (from X = 0 to X ∼ 20 Å in Figure ) was used to extract the bulk propane density. The bulk pressure for each system was then estimated from the bulk pure propane density in gas phase by using the Peng–Robinson equation of state.[41] A schematic of the process is provided in Figure . The estimated bulk pressures for five C3H8–H2O mixtures of different compositions are given in Table .

Figure 2

Detail of a simulation snapshot revealing the distribution of propane away from the solid substrate (a), and propane molecular density profiles along the X direction of the simulation box, in the region outside the cylindrical pore (b). In this representation, the solid substrate is located at x > ∼ 46 Å. Note that the density profiles show accumulation of propane near the substrate, as visualized in the snapshot. The vertical line in panel a identifies the bulk region within which the propane density is constant. The white lines on the right of panel a help identify the silica substrate, which yields a cylindrical pore filled with water (red and white for O and H atoms, respectively).

Table 2

Bulk Propane Pressures Estimated for Five Propane–Water Systems Simulated within the System A Model at T = 300 K, and the Composition of the Corresponding Systems Simulated with System B Models

number	System A composition	estimated bulk pressure (MPa)	System B composition
1	76 C3H8–2141 H2O	0.60 ± 0.05	11 C3H8–387 H2O
2	176 C3H8–2141 H2O	1.5 ± 0.1	15 C3H8–379 H2O
3	221 C3H8–2141 H2O	1.9 ± 0.1	17 C3H8–383 H2O
4	281 C3H8–2141 H2O	2.30 ± 0.02	21 C3H8–373 H2O
5	350 C3H8–2141 H2O	2.90 ± 0.05	22 C3H8– 371 H2O

System B

To quantify the properties of confined fluids, the resulting configurations obtained from the simulations conducted using the System A setups were modified by removing the region outside of the pores (along the X direction), and by rendering the cylindrical pores effectively infinite. For each system, the number of propane and water molecules inserted inside the pore corresponded to those obtained from the System A simulations. In Table we report the number of fluid molecules confined within the corresponding System B pores.

To further quantify the impact of water on the diffusion of confined propane, we conducted additional simulations in which (a) 45 propane molecules were confined within the cylindrical pore (no H2O present, corresponding to bulk pressure of ∼0.6 MPa), and (b) starting from system B5 (22 C3H8 and 371 H2O molecules) we systematically reduced the amount of water molecules (three simulations with 321, 271, and 221 water molecules, respectively, while maintaining 22 C3H8 molecules confined within cylindrical pores).

Force Fields

For all simulations in which fluid molecules were at contact with the solid substrate, the CLAYFF[42] force field was implemented to simulate the silica substrates. CLAYFF is a general force field suitable for fluid-clay and other clay-related systems. In these simulations, the silica frame was kept rigid by freezing the position of the bulk O and Si atoms while only the surface H atoms of silanol groups were allowed to move. Because the simple point charge/extend (SPC/E) model[43] provides acceptable estimates for the structure, the internal energy, the density and the diffusivity for water at ambient conditions,[44] it was selected to describe water. Building on prior simulation results from our group,[6] the TraPPE-UA force field[45] was employed to model propane. Dispersive forces were described implementing the 12–6 Lennard–Jones potential and electrostatic forces were taken into account for nonbonded interactions. The distance cutoff of interatomic interactions for all simulations was fixed at 14 Å. The particle mesh Ewald (PME) method[46] was used to treat the long-range electrostatic interactions. The Lorentz–Berthelot mixing rule[47] is applied to determine the Lennard-Jones parameters for unlike interactions.

Algorithms

All simulations were carried out using the Groningen Machine for Chemical Simulations (GROMACS) simulation package, version 5.1.2.[48,49] The leapfrog algorithm[50] was used to integrate the equations of motion. The temperature of the silica matrix was kept constant at 300 K using the Nosé–Hoover thermostat. The temperature of the confined fluids was also maintained constant at 300 K using the Nosé–Hoover thermostat. It was found that decoupling the two thermostats prevents unrealistic distributions of the kinetic energy between solid and fluid.[51] Both thermostats had a fixed temperature-damping factor of 100 fs. Following the minimization of the energy for the initial configuration, both systems A and B were simulated using molecular dynamics.

For Systems A, we found that equilibration was achieved after ∼60–100 ns of simulation time, depending on loading. Equilibration was considered achieved when the propane densities oscillated around constant values, and both system temperature and energy fluctuations remained within 10% of their respective average values. The density profiles of propane and the adsorption isotherm obtained from the last 2 ns of the simulations are presented below as results.

For Systems B, all systems were simulated for 80 ns with a time step of 1 fs. After 78 ns of equilibration time, trajectories from the last 2 ns of simulations were used to analyze transport and structural properties of the fluids. Equilibration was considered achieved when the density profiles of propane and water within the pore and the total energy of the systems converged within the criteria discussed above. The self-diffusion coefficients were calculated from the mean square displacement by implementing the Einstein equation:[52]In eq , r(t) and r(t′) are the positions of particle i at time t and at the time origin t′, respectively, and d is the number of degrees of freedom. The diffusivity of fluids within the pore is considered as a one-dimensional translation along the X direction, because the cylindrical shape of the pores constrains the movement of molecules over long distances along the Y and Z directions.

Simulation Results

Adsorption Isotherms

In Figure , we report the simulated amount of propane adsorbed in the cylindrical pores filled with water at 300 K and bulk pressures in the range ∼0.6 to 2.9 MPa. For these simulations, we used the approach described as System A in the Simulation Methods and Algorithms section. The amount of propane adsorbed in the pore increases as the pressure increases. For comparison, the experimental mole fraction solubility of propane in bulk water at 298 K and 1 bar is 2.732 × 10–5,[53] which corresponds to ∼0.067 g of propane per kg of water. For completeness, we point out that Ferguson et al.[54] used molecular dynamics to study the solubility of linear alkanes in water in the bulk. They found excellent agreement between simulated and experimental solubility for short alkanes. Ferguson et al. implemented the SPC/E force field to simulate water, and the TraPPE force field to simulate alkanes. Based on these observations, the results in Figure suggest that the amount of propane in confined water is much larger than that which could be expected based on bulk solubility data. However, analysis of the simulation snapshots, discussed below, show that in the case simulated here propane is not solubilized in confined water.

Figure 3

Simulated amount of propane adsorbed in the cylindrical pore filled with water as a function of the bulk pressure. The simulations were conducted at 300 K. The amount of propane adsorbed is expressed as grams of propane per kilogram of water inside the pore.

Structural Properties: Atomic Density Profiles

Atomic density profiles as a function of distance from the pore axis were calculated for all confined fluid molecules. These simulations were conducted using the approach described as System B in the Simulation Methods and Algorithms section, but they were initiated from data obtained from the final configurations derived from simulations for Systems A. All simulation results described in the remainder of the manuscript were obtained using the System B set up. Representative density profiles are shown in Figure , where the density profiles of silanol group atoms are used to locate the pore surface, and different lines represent density distributions of various compounds. We report the atomic density profiles of O and H atoms of water and the molecular density profile of propane, as obtained from the distribution of the CH2 group of confined propane. The density profiles of confined water do not change significantly as the amount of propane increases, hence only one data set is shown for O and H atomic density profiles.

Figure 4

(a) Axial view of a representative simulation snapshot illustrating the distribution of fluid molecules confined in the cylindrical silica pore. (b) Radial density profiles calculated for molecules within the cylindrical silica pore. The reference 0 is the central axis of the pore. For water we report the density profiles of both oxygen and hydrogen atoms. These density profiles do not change considerably with pressure, as the amount of water in these systems does not change significantly. For propane, the density profile is obtained from the position of the ethyl pseudoatoms.

The radial density profiles for the O atom of water, shown in Figure , are characterized by one pronounced peak at ∼7.4 Å. This peak indicates the formation of one layer of water in proximity of the solid substrate. The formation of this hydration layer seems to be in agreement with experimental neutron scattering data obtained for water confined in MCM-41.[55,56] The density profile for H atoms of water provides further information. In particular, two peaks are observed around the position of the pronounced O density peak: one at ∼7.1 Å and one at ∼8.4 Å. Because the H peak at ∼8.4 Å is approximately identical in intensity as the O peak at ∼7.4 Å, the results suggest that the water molecules in the first hydration layer lay with one of their OH vectors toward the surface. The second H peak, located at ∼7.1 Å has nearly twice the intensity than the O peak at ∼7.4 Å, suggesting that this peak is due to H atoms of water molecules in the first layer plus H atoms of water molecules whose oxygen atoms are in the second hydration layer. The peaks in both O and H density profiles at distances larger than ∼8.4 Å represent water molecules that are found within the atomically rough surface. The atomic density profiles of confined water show that water does not homogeneously occupy the pore volume, as water seems to be depleted from the pore center. Not far from the pore center, the water density approaches the bulk value (0.033 molecules/Å3). However, given the small pore and the features already discussed, water in this region probably does not behave as bulk water. Few water molecules are found at r > 8 Å due to the surface roughness of the substrate. These molecules (about 22) are found to penetrate small cavities present within the silica matrix, and are excluded from the calculations of the dynamical properties of confined water presented below. The water atomic density profiles do not change significantly as the bulk pressure increases, since the number of water molecules within the pore is not considerably different.

On the other hand, increasing the bulk pressure by adding propane molecules to the simulations of System A yields an increase in the molecular density of propane within the pores. Note that propane density peaks, irrespective of the simulated pressure, are always close to the pore center, implying that propane molecules preferentially accumulate in this region, and remain excluded from the hydration layer.

The preferential orientations of confined molecules with respect to the radial direction are reported in Figure . For propane, we quantified the angle θp formed between the CH3–CH3 vector identified by each propane molecule, and the vector pointing from the pore surface to the pore center. For water, we quantified the angle θw formed between the dipole moment vector of a water molecule and the vector pointing from the pore surface to the pore center. Schematics of these angles are illustrated in Figure a. For propane molecules, when the angle θp is 0° or 180°, the CH3–CH3 vector is perpendicular to the pore surface, whereas, when θp is 90°, the propane lays parallel to the surface. The interpretation for the results on the orientation of the dipole moment of water molecules is analogous. Our results in Figure b show that propane molecules preferentially orient parallel to the pore surface, irrespective of the distance from the surface. Conversely, water molecules show different orientations along the radial distance. In the central region of the pore, the dipole moment of H2O molecules yields a preferential angle of ∼90° with respect to the radial direction, while the preferential angle shows a strong dependence on radial position near the pore surface. The results obtained for θw in the interfacial region suggest that one OH group of the water molecules points generally toward the surface. This result is consistent with the hydrogen-up orientation of water molecules observed for the atomic density profiles of water and discussed above (see Figure ).

Figure 5

Scheme representing the orientation angles θw and θP as calculated for confined water and propane molecules, respectively (a). Orientation of confined fluids as a function of their distances from the central axis of the pore (b) for a system composed of 11 propane molecules and 387 water molecules.

To quantify the structure formed by water molecules within the cylindrical pore, we computed the number of hydrogen bonds (HBs) that one water molecule forms with other water molecules or with the silica surface. The results are displayed in Figure in the form of the density of HBs (panel a), or of the average number of HBs per water molecule (panel b), wherein both cases are a function of the radial distance. The geometric criterion proposed by Martí[57] was implemented to determine when a HB is formed. According to this criterion, one HB is formed when the distance between the acceptor oxygen and the donor hydrogen is less than 2.4 Å and the H–O···O angle between the atoms involved in the HB is lower than 30°. It was found that changing the amount of propane has little effect on these results. Only one data set is shown for clarity. For reference, when this criterion is used to assess the number of HBs per water molecules in bulk liquid water, it is found that one water molecule forms on average 3.4 HBs.[57] The results for the HB density profile (panel a in Figure ) show a pronounced peak located at ∼5.5 Å. The water molecules in this position (see Figure , panel b) adopt an angle θw of ∼113°. This implies that water–water HBs dominate within this region, where water molecules are too far from the surface for water molecules to form HBs with the surface silanol groups. This is confirmed by the data shown in panel b of Figure , which further show that water molecules in this region form an average of ∼1.5 HBs per molecule. By contrast, the water molecules in the interfacial region (radial position of ∼8 Å) preferentially form HBs with silanol groups, in some cases yielding ∼2.8 water-silanol HBs per water molecule. Comparing the data on HB density profiles and those for the preferential orientation of water molecules, we conclude that the oxygen atoms of interfacial water molecules tend to serve as acceptors for HBs formed with the hydrogen atoms provided by the surface silanol groups, which act as donors.

Figure 6

(a) Density profile of hydrogen bonds. (b) Average number of hydrogen bonds per water molecule. The data are plotted as a function of the radial distance from the central axis of the pore. Hydrogen bonds are distinguished as forming either among water molecules or between water molecules and surface silanol groups. The results presented in this figure were obtained for a system composed of 15 propane and 379 water molecules.

Dynamical Properties 1: Translational Dynamics of Confined Propane

The self-diffusion coefficients for propane are summarized in Table . For reference, the self-diffusion coefficient for propane has been simulated by Feng et al.[58] using the OPLS-UA model. At 294 K and 25 MPa, the simulated value was ∼9.08 × 10–9 m2/s, which compares well with the corresponding experimental value of 9.095 m2/s measured by NMR spin echo by Greiner-Schmid et al.[59] The results in Table show a drop in propane self-diffusion coefficient as pressure increases. The self-diffusion coefficient for confined pure propane was also computed for comparison, and it is reported in Table as well. The corresponding result, which is comparable to the self-diffusion coefficient for bulk propane, yields a considerably larger 1D self-diffusion coefficient compared to that obtained in the presence of water at the same bulk pressure (by a factor of ∼23). These results confirm that the presence of water strongly impedes the transport of propane across the hydrated pores, in qualitative agreement with prior simulations,[60] and also with expectations. A similar damping of propane mobility due to presence of D2O in MCM-41-S has also been observed in quasielastic neutron scattering experiments.[31] Unfortunately, only a qualitative comparison can be made, because the experiments were conducted at lower temperatures than the simulations presented here.

Table 3

One-Dimensional (1D) Self-Diffusion Coefficient Estimated for Propane Confined in Silica Pores at Different Pressuresa

	D (10–10 m2/s)
	0.6 MPa	1.5 MPa	1.9 MPa	2.3 MPa	2.9 MPa
confined propane with H2O	4.1 ± 0.2	3.0 ± 0.2	2.8 ± 0.1	2.4 ± 0.1	1.8 ± 0.1
pure confined propane	94.5 ± 0.8

Errors are estimated as one standard deviation from the average.

Dynamical Properties 2: Residence Times

To quantify the dynamic properties of water in the systems considered, we computed the residence autocorrelation function CR(t). This quantity allows us to estimate how long a water molecule, found at a specific location away from the surface, remains in that position. The algorithms are described elsewhere.[61−63] In general, the faster CR(t) decays from 1 to 0 as time progresses, the faster molecules leave the layers considered for these analysis. The analysis is focused on water molecules that reside within two water layers of interest, denoted as layer I (radial position >6.8 Å) and layer II (5.0 Å < radial position <6.5 Å). In each case, water found within an annular region of thickness 1.5 Å was considered. The results are reported in Figure . It was found that changing the amount of propane has little effect on these results. Only one data set is shown for clarity. We also report in Figure analogous autocorrelation function results obtained for water molecules on a flat crystalline silica substrate as reported by Ho et al.[64] Our results show that water molecules at contact with the pore surface stay within the hydration layer (layer I) longer than water molecules found in layer II stay in that region. This is probably due to the preferential interactions between water and silica surface, perhaps via hydrogen bonds. Further, we observe that water molecules reside for longer time in the hydration layer within the cylindrical pore than on the crystalline flat substrate. These results suggest that the shape and the surface properties of the support strongly affect the dynamical properties of water in the hydration layers.

Figure 7

Residence autocorrelation functions CR(t) for oxygen atoms in the first and the second hydration layer of the hydrated silica pores. Solid lines represent results obtained for water confined in the cylindrical amorphous pore considered here, while dashed lines represent results obtained for water on a crystalline flat silica substrate with surface density of 4.54 OH/nm2 reported by Ho et al.[64] Results indicated as LCO-1 and LCO-2 are for water molecules found within a hydration layer centered at 0.95 and 2.45 Å from the flat surface, respectively. The results for Layer I and Layer II were obtained for a system composed of 17 propane and 383 water molecules.

Dynamical Properties 3: Rotational Dynamics

The rotational dynamics of molecules can be quantified by calculating vector–vector autocorrelation functions. These calculations reveal changes over time in the orientation of selected molecules. Following previous reports,[61−63] we computed reorientation autocorrelation functions defined asIn eq , v(0) is either the dipole moment, the hydrogen–hydrogen vector (HH), or the CH3–CH3 vector of molecule i at time t = 0; v(t) is the same quantity for molecule i, at time t.

The results in Figure a suggest that the CH3–CH3 vector autocorrelation functions decay more slowly when increasing bulk pressure, which is attributed to steric hindrance. This observation is in agreement with the results discussed for the diffusivity of propane molecules, which decreases as the bulk pressure increases. We estimate the time required by the CH3–CH3 vector autocorrelation function to decay from 1 to 1/e, which was found to be 1.82, 2.62, 3.50, 4.00, 4.62 ps, respectively, as P increases from ∼0.6 to ∼2.9 MPa. In Figure b,c, we report relevant autocorrelation functions estimated for confined water. It was found that these results do not depend strongly on the amount of propane present (the exception being the dipole–dipole autocorrelation function in which case some difference was observed). For brevity, only the results from one system are shown. Both dipole–dipole and HH vector autocorrelation functions in layer I decay more slowly than they do in layer II. These results indicate that those water molecules found in layer I, which are highly associated with surface hydroxyl groups, have a slow rotation, while water molecules in layer II have higher reorientation freedom. These results are qualitatively consistent with those reported by Milischuk et al.[20] Clearly, solid–water interactions strongly impact the rotation of water molecules. We further note that the reorientation autocorrelation function of the HH vector decays faster than that of the dipole moment vector, suggesting that the rotation of water in the layers considered is anisotropic. Finally, the reorientation dynamics of water molecules confined within the cylindrical pore considered here are found to be much slower than that observed for hydration water on a crystalline flat silica substrate.[12] This is surprising, since in the crystalline flat substrate the OH bonds in the substrate were treated as rigid, while they are allowed to vibrate in the present work. This suggests that the cylindrical morphology of the surface effectively enhances the strength of the preferential interactions between water molecules and the solid substrate. This is consistent with the high number of HBs formed per water molecule at the interface with the solid substrate (see Figure ).

Figure 8

CH3–CH3 vector autocorrelation function for confined propane molecules at different bulk pressures (a), dipole moment (b), and hydrogen–hydrogen vector (c) autocorrelation functions for water molecules in layer I and layer II within the hydrated silica pores. Solid lines represent results obtained for water in the cylindrical pores considered here for a system composed of 15 propane and 379 water molecules. Dashed lines represent results obtained for water on a flat crystalline silica substrate with hydroxyl surface density of 4.54 OH groups per nm2 reported by Argyris et al.[12] Results indicated as BO-1 and BO-2 are for water molecules found within a hydration layer centered at 2.15 and 3.05 Å from the flat surface, respectively.

Hydrogen Bond Network: Dynamical Properties

We assessed the average lifetime of HBs as a function of the position within the cylindrical pore by calculating the HB autocorrelation function, defined as[65]In eq the variable h(0) equals 1 when a HB is found at time t = 0. If the tagged HB remains formed as the time t progresses, then h(t) remains equal to 1. The quantity h(t) switches to 0 when the HB breaks, and it remains 0 afterward. The results presented in Figure are averages of 5 calculations. It was found that changing the amount of propane has little effect on the results obtained for confined water. Only one data set is shown for clarity. They show that the HB autocorrelation function for water molecules in layer I decays more slowly at short time intervals than that obtained for water molecules in layer II. This indicates that the water–water HBs near the pore surface remain intact longer compared to those further from the surface. Similarly, we found a significantly slow decay for the HBs established between water and surface silanol groups (within layer I). This result suggests that water molecules in layer I prefer to form stable HBs with the pore surface rather than with other water molecules. This is expected when considering that interfacial water molecules show slower CR(t) and rotational diffusion when compared to those further from the pore surface, as discussed above. The results for HB–HB autocorrelation function of the HBs formed between water molecules near a crystalline flat substrate reported previously by Argyris el at.[61] are also shown in Figure . Consistent with the dynamical results obtained so far, our data suggest that water–water HBs last longer when water is near the cylindrical pore surface than near the flat substrate. Because the density of −OH groups is much larger on the crystalline substrate considered by Argyris et al. than it is for the cylindrical pore considered here, it seems that the cylindrical pore geometry has a stronger effect on these autocorrelation functions than the density of −OH groups.

Figure 9

Hydrogen bond–hydrogen bond autocorrelation function for water molecules found within layer I and layer II in the cylindrical pores considered here (solid lines). These results were obtained for a system composed of 21 propane and 373 water molecules. The dashed orange line corresponds to hydrogen bond–hydrogen bond autocorrelation function of hydrogen bonds formed between water molecules on the partially hydroxylated slab pore surface with a total surface density of 6.8 −OH/nm2 reported by Argyris el at.[61]

Effect of H2O Loading on Transport of Confined Propane

Both dynamical and structural features used to quantify the properties of confined fluids suggest the formation of a stable, packed hydration layer (radial distance >6.8 Å), with a more diffuse and perhaps sparsely filled region of water near the center of the pore. Our results suggest that propane molecules accumulate in this region, near the center of the pore, where water molecules can be depleted. The results for the transport properties of confined propane (see Table ) suggest that the translational diffusion is much slower in the hydrated pores than in pores filled only by propane. In Figure we compare the CH3–CH3 autocorrelation function for propane calculated when propane is the only fluid in the pores, and when water is also present. The results suggest a significant difference between the rotational diffusion of confined propane. The time required for the CH3–CH3 vector autocorrelation function to decay from 1 to 1/e in the hydrated pore is ∼3 times longer than that obtained in the dry pore. It is perhaps interesting to point out that this delay in the rotation of propane molecules is not due to its interactions with the solid matrix, but rather to its interactions with water molecules found within the pore.

Figure 10

CH3–CH3 vector autocorrelation functions for propane molecules. When only propane is present, 45 molecules are considered, while when the pore is hydrated we considered 11 propane molecules and 387 water molecules.

The results presented so far suggest that the transport properties of confined propane are controlled by the large number of water molecules present in the pore. To test this hypothesis, we assessed the effect of water loading. We calculated atomic density profiles and self-diffusion coefficients for both water and propane for systems in which the propane loading was maintained constant and the amount of water was reduced. The atomic density profiles are shown in Figure . In each simulation, water density profiles show that water always accumulates near the silica substrate, yielding a hydration layer consistent with layer I in the analysis above. The peak density for water within the hydration layer decreases as the water content is reduced. In other words, as more water fills the pore, the first hydration layer builds up. At low hydration levels, almost no water molecules are found near the pore center. In all cases, the density profiles for propane accumulate near the pore center, suggesting that reducing H2O loading does not promote the adsorption of propane on the silica surface at least for the systems considered here. This was expected, given the low amount of propane, the large amount of water, and the preferential interactions between the silanol surface groups and the water molecules.

Figure 11

Radial density profiles of molecules confined within the amorphous cylindrical silica pore. Systems 1, 2, and 3 contain 321, 271, and 221 water molecules, respectively. In all cases, the number of propane molecules is kept constant at 22.

We report in Table the one-dimensional self-diffusion coefficients for the systems considered in Figure . We computed the self-diffusion coefficients for only those water molecules found in the middle region of the pore (radial distance <6 Å), since the water molecules in the hydration layer form stable HBs with the surface.

Table 4

1D Self-Diffusion Coefficient Estimated for Confined Fluids in Silica Pores for Four Systems with Different Water Loading

		self-diffusion coefficient (10–10 m2/s)
system	system composition	C3H8	H2O
B5	22 C3H8–371 H2O	1.8 ± 0.1	1.55 ± 0.05
1	22 C3H8–321 H2O	8.7 ± 0.1	2.0 ± 0.1
2	22 C3H8–271 H2O	27.9 ± 0.5	1.48 ± 0.13
3	22 C3H8–221 H2O	64.2 ± 0.8	0.6 ± 0.1

The self-diffusion coefficient of confined propane was found to decrease as the water loading increases; conversely, the self-diffusion coefficient of confined water was found to increase, reach a maximum, and then decrease. Analysis of sequences of simulation snapshots confirms that the reduced propane diffusion is due to the formation of “molecular bridges” formed by water molecules across the silica pore, which hinder the free flow of propane transport. In Figure we report a snapshot for such a molecular bridge.

Figure 12

Simulation snapshots representing the flow patterns of water molecules across the pores. Panel a shows bridges of water molecules formed within the pore filled 22 propane molecules and 371 water molecules. Panel b shows the dissolution of the molecular bridges within the pore filled with 22 propane and 271 water molecules.

For the systems considered in Table , the molecular bridges are not present when 321 or fewer water molecules are present, but they form when 371 water molecules are simulated. Further, a decrease in self-diffusion coefficient of water as water loading decreases below 321 corroborates the formation of stable hydration layer near the pore surface as discussed above. The significant decrease in self-diffusion coefficient when water loading increases from 321 to 371 is therefore mostly due to the formation of water bridges, which hinder the translational diffusion of water molecules found near the pore center. In fact, these molecules are now hydrogen-bonded to water molecules all around the cylindrical pore.

Conclusions

A series of atomistic molecular dynamics simulations were performed for systems composed of water and propane confined within cylindrical silica pores. The simulations were designed to study the effect of bulk pressure and water loading on the mobility of confined propane. The pressure was controlled by changing the amount of propane. The pore was a model amorphous cylindrical silica pore designed to resemble the pores in MCM-41 materials. All simulations were conducted at 300 K. The simulation results are quantified by analysis of the composition of the confined fluid systems, molecular density profiles in the radial direction, preferential orientations of the fluid molecules with respect to the pore, and both rotational and translational dynamical properties. Our results reveal that propane accumulates near the pore center, where water can be depleted. Conversely, water molecules tend to form hydrogen bonds with the silanol groups on the pore surface. We found that the self-diffusion coefficient of confined propane decreases as bulk pressure or water loading increase. The significant effect of water on the diffusion of confined propane is due to the formation of water bridges that span the pore volume, thus hindering propane transport. This observation provides molecular-level interpretation for recently reported experimental findings regarding mixtures containing propane and D2O confined in MCM-41-S. The qualitative agreement between simulations and experiments contribute to our understanding of transport of hydrocarbons in subsurface environments.