Structure and Dynamics of Interfacial Water on Muscovite Surface under Different Temperature Conditions (298 K to 673 K): Molecular Dynamics Investigation

Abstract

We performed molecular dynamics (MD) simulations to study structure, stability, and dynamics of the water adsorption layer on muscovite mica at several temperatures (from 298 K to 673 K) and pressures (0.1 MPa, 10 MPa, and 50 MPa). We studied the structure of the adsorption layers with three characteristic peaks of density and orientation of H₂O molecules in one-dimensional and two-dimensional profiles. The results show that the water adsorption layers become less structured and more mobile as the temperature increases. We also found the first and the second layers are less diffusive than the third one, and the difference of diffusivity gets unclear as the temperature increases. Finally, we discuss implications to hydration forces and wettability, which are significant interfacial properties of the multiphase fluids system such as water/gas/mineral systems, from the viewpoint of water adsorption film with nanometer thickness.

1. Introduction

In subsurface porous media with high specific surface area, water–mineral interfacial property is a key factor for mass transport (e.g., [1,2,3]). Better understanding of this is beneficial not only for revealing mechanisms of natural phenomena but also for improving the effectiveness of subsurface engineering processes such as CO₂ geological storage (CGS), geothermal energy production, and oil/gas recovery.

For example, wettability, which is affected by fluid–mineral interactions, is significant for the seal integrity of CGS (e.g., [4,5,6]), and the recovery ratio of oil production [7,8]. Permeability, which is also affected by water–mineral interaction (e.g., [9]), is a critical factor for recovery of geothermal energy [10,11,12]. Moreover, the water–mineral interaction reduces the amount of CH₄ adsorption on shale, which largely affects the potential of shale gas production [13,14]. Furthermore, the water adsorption film in water/gas unsaturated system is an important media for mass transport [2,15,16].

Interfacial water on the mineral surface region, of which thickness is from micro- to nanometer, shows different structure and properties from that in bulk phase (e.g., [17,18]). These differences have been studied by using experimental, theoretical, and numerical approaches. For example, X-ray reflectance measurements and surface force measurements using atomic force microscopes have clarified nanometer-order structures [19,20,21,22]. It has been reported that the distribution of water molecules above the muscovite surface exhibits three distinct peaks [19,20,21,22]. Results obtained from molecular dynamics (MD) simulation are in good agreement with a measured water adsorption structure, and more and more research using MD simulation has been conducted (e.g., [23,24,25,26,27,28,29,30,31]).

Water/muscovite interfacial properties have been studied by experiments and MD simulations; however, results of the previous studies focused on the ambient condition. Although the temperature change is considered to affect the properties of water adsorption layer on the mineral surface [32,33], no quantitative investigation on its structure and dynamics has been conducted under high temperature conditions such as 673 K. Considering the geothermal gradient, it is crucial to study the interfacial properties in a wide range of temperature conditions. For example, 10 MPa/313 K is one of the typical reservoir conditions for CGS [34]. A much higher pressure–temperature condition such as 50 MPa/673 K needs to be considered for supercritical geothermal resources recovery [35].

Therefore, in this study, MD simulation of the water/muscovite interfacial system is performed at several temperatures (from 298 K to 673 K) and pressures (0.1 MPa, 10 MPa, and 50 MPa) focusing on the changes in the structure and the dynamical properties of the water adsorption layer. This paper is organized as follows. In Section 2, calculation details in this study are explained. In Section 3, we show the structural and dynamical properties of the interfacial water on muscovite. We discuss implications of the results to hydration forces and wettability in Section 4. Finally, the paper is summarized in Section 5.

2. Materials and Methods

We performed MD simulation of a water/muscovite system under several pressure and temperature conditions shown in Table 1. Water becomes supercritical at 50 MPa/673 K, while it is liquid state under the other conditions [36].

In this study, we did not focus on the chemical reactions on the muscovite surface. Although the melting point of muscovite is higher than 1500 K, which is much higher than the temperatures considered here, dissolution of muscovite in water and its temperature dependence has been reported (e.g., [37]). In this study, MD simulations were performed to predict the structure and dynamics under high temperature and pressure conditions, validity of which was confirmed by referring to the previous MD results under ambient conditions.

First, the calculation system was stabilized energetically. Then, physical properties, such as density and energy, were evaluated as ensemble averages. In the preliminary run, the temperature and pressure of the system were maintained by the Berendsen thermostat and barostat [38]. After equilibration, NPT MD calculations were performed for 100 ns (1.0 × 10⁸ steps) using the Nosé–Hoover chain thermostat [39,40] (the time constant ${\tau }_T$ = 0.1 ps) and Parrinello–Rahman barostat [41] (the time constant ${\tau }_p$ = 5.0 ps). The physical properties were estimated using the last 50 ns data (5.0 × 10⁷ steps) of the NPT run in equilibrium. The interatomic interaction $U_{ij}$ was calculated based on the Lennard Jones (LJ) potential and the Coulomb potential, as shown below:

(1)$U_{ij}=4{\epsilon }_{ij}\left[{\left(\frac{{\sigma }_{ij}}{r_{ij}}\right)}^{12}-{\left(\frac{{\sigma }_{ij}}{r_{ij}}\right)}^6\right]+\frac{1}{4\pi {\epsilon }_0}\frac{q_i{q}_j}{r_{ij}}$

where ${\epsilon }_{ij}$ and ${\sigma }_{ij}$ are the energetic and size parameters for the LJ potential, respectively: $q_i$ and $q_i$ are the charge of particle $i$ and $j$: $r_{ij}$ is the distance between particle $i$ and $j$: ${\epsilon }_0$ is the dielectric permittivity of vacuum (8.85419 × 10⁻¹² F/m).

The parameters of non-bonded interactions in the Equation (1) for H₂O and muscovite were determined after the SPC/E model [42] and the CLAYFF model [43], respectively. Terms of bonded interactions are incorporated in these models for describing the energy of stretch of the hydroxyl bond and angle bend of hydroxyl groups. Since parameters of the CLAYFF are developed for reproducing observed physical properties of clay minerals including interface with water, it has been widely used for MD simulations of water/clay systems (e.g., [22,23,24,27,28,29,30,31]). Since its parameters were determined based on the SPC-type water model, we chose the SPC/E models as the force field of H₂O molecules for consistency with the CLAYFF. The Lorentz–Berthelot rule was used for ${\epsilon }_{ij}$ and ${\sigma }_{ij}$ of unlike atoms. The cutoff radius was 1.1 nm for both the LJ potential and the short-range Coulomb potential. The long-range Coulomb potential was calculated by the particle mesh Ewald (PME) method. GROMACS 2019 [44] was used for MD simulations. VMD (version 1.9.3) [45] was used to visualize the snapshots.

Calculations were performed using the system shown in Figure 1. H₂O phases are in contact with (0 0 1) plane of the muscovite surface, which was set to be orthogonal to the Z-axis of the simulation box. The area marked by the blue line in Figure 1 is the original calculation system, and the area outside the frame is the image cell duplicated by the periodic boundary conditions. This results in a system of slit pores in which the fluid phase is sandwiched between the mineral phases. Muscovite is a 2:1 layered clay mineral, and its chemical composition is represented by KAl₂AlSi₃O₁₀(OH)₂. The structure of the muscovite model is based on the X-ray measurement [46]. Muscovite has a stable cleavage surface on (0 0 1) plane (Figure 2a) and is commonly used for surface force measurement experiments. The tetrahedral layer and the octahedral layer are included at a ratio of 2:1 (Figure 2b,c). The mineral surface is negatively charged by isomorphous substitution of Si for Al at the tetrahedral layer at a ratio of 3:1 (red dotted line in Figure 2b,c). Cations are adsorbed on the charged surface and a K⁺-muscovite model was adopted in this study. Hereinafter, adsorption sites of K⁺ are called K⁺-site, and the other sites are called o-site as seen in Figure 2d. The same arrangement as previous studies [24,25,26,27,28] was used for the isomorphic substitution.

The unit cell size of muscovite is (ΔX, ΔY, ΔZ) = (0.5199 nm, 0.903 nm, 2.011 nm), and the unit cell was duplicated 6, 3, and 1 times in the X, Y, and Z directions, respectively. Therefore, the dimension of the muscovite phase in the model was (3.119 nm, 2.709 nm, 2.011 nm). In a closed slit pore system, dynamic behaviors and hydration forces in confined water can be different from those in an open system, as reported by Calero and Franzese (2020) [47]. In this study, the water phase consisting of 1000 H₂O molecules (for the systems at 0.1 MPa/523 K and 0.1 MPa/573 K) or 1400 H₂O molecules (for the rest of the systems) is in contact with the muscovite surface. This resulted in a distance between the surfaces large enough to evaluate density and diffusivity relative to the bulk phase, under the pressure and temperature conditions in this study. The distance at each system is shown in Table S1 in Supplementary Materials.

3. Results and Discussion

3.1. Structural Properties

3.1.1. One-Dimensional (1D) Density Profile

1D profiles of water molecules as a function of distance from the muscovite surface are shown in Figure 3. Graphs normalized by the bulk density for each water molecule, oxygen atom, hydrogen atom, and K⁺ are shown in Figures S1–S5 in Supplementary Materials. At 0.1 MPa and 298 K, the distribution has three characteristic peaks near the muscovite surface. This is in good agreement with the result of measurement data [19,20,21,22,23] and existing MD simulations [22,23,24,25,26,27]. As the temperature increases, density values decrease, both near the muscovite surface and in the bulk region. The height of the first peak is not significantly affected by increasing temperature. However, the height of the second peak decreases with increasing temperature. Moreover, the third peak gradually becomes unclear with temperature.

A similar trend is seen not only at 0.1 MPa but also at higher pressures of 10 MPa and 50 MPa. At 50 MPa and 673 K, at which water is in the supercritical state, the height of the adsorption layer relative to the bulk density is higher than any other conditions, as shown in Figures S1–S3 in Supplementary Materials. Moreover, the distinction of each individual peak becomes unclear. It can be said that the adsorption region tends to lose its characteristic structure as the temperature increases.

3.1.2. Two-Dimensional (2D) Density Map

Here, a detailed analysis was performed to gain more insights into the three characteristic peaks seen in Figure 3. According to the distance from the muscovite surface $d$ nm, the first, second, and third layers are defined as the region having the width of 0.1 nm <$d$ < 0.2 nm, 0.2 nm < $d$ < 0.5 nm, and 0.5 nm < $d$ < 0.7 nm, respectively. For these three regions, the density maps on the XY plane were calculated as shown in Figure 4. To clearly show the density of adsorbed water on the (0 0 1) surface of muscovite, the lattice is also overlaid on the density maps. For comparison, the density values on the maps are normalized by bulk density as in Figures S2 and S3 in Supplementary Materials. Two-dimensional density profiles of (a) the first layer, (b) the second layer, and (c) the third layer at four different conditions (0.1 MPa/298 K, 50 MPa/298 K, 50 MPa/473 K, and 50 MPa/673 K) are shown in each panel of Figure 4, respectively. The detailed characteristics of each layer are stated below.

• The first layer

At the center of the hexagonal lattice, high-density areas of oxygen are seen. This is consistent with the 2D density map calculated by the previous MD study [29]. It is found that O and H of H₂O can reside only at o-site. Figures S4 and S5 in Supplementary Materials show that the height of K⁺ adsorption is about 0.15 nm from the muscovite surface, which is in good agreement with previous measurement results (0.167 nm–0.177 nm by Schlegel et al. (2006) [48]). In addition, as shown in Figure S7 in Supplementary Materials, K⁺ resides at the center of the hexagon, which is in good agreement with a previous study by Ricci et al. (2014) [49]. Since the adsorbed K⁺ on the surface is found in the first layer, H₂O molecules adsorb at different sites other than the K⁺-site. This is also seen in the snapshot in Figure 5b. As the temperature increases, high-density areas of hydrogen atoms become wider. This indicates that H₂O molecules at the adsorption site vibrate more intensively at high temperatures than at lower temperatures, although the height of the first peak in Figure 3 is almost unchanged.

• The second layer

On K⁺-site, the water density is lower than in the other regions. In order to understand this, the density profiles of O and H on K⁺-site and o-site are shown in Figure 6. For calculating these local density profiles, K⁺-site is defined as the region within 0.138 nm (ionic radius of K⁺) from the position K⁺ is adsorbed on the muscovite surface (see Figure 2 and Figure S8 in Supplementary Materials), while the other sites are defined as o-site. In Figure 6, the profiles on o-site are similar to those in Figure 3. In contrast, on K⁺-site, neither O nor H exists around K⁺ ($d$ < 0.3 nm; the first and second layers), which is quite different from the profiles on o-site. The low-density areas on K⁺-site on the 2D map in the second layer in Figure 4b can be interpreted from these local density profiles in Figure 6.

• The third layer

High-density areas of O and H of H₂O are seen even at K⁺-site, although the density profiles become more homogeneous than those in the first and second layers. As the temperature increases, the density is more homogeneously distributed. This is consistent with the trend in Figure 3 that the third peak becomes less sharp as the temperature is increased.

3.1.3. 1D Profile of the Orientation Angle

Figure 3 suggests that there exist three adsorption layers of water on the muscovite surface. Previous studies reported that the formation of distinct water layers is related to the polarity-based orientation of the water molecules [22]. To discuss the effect of the temperature on the structural change, the orientation angle was calculated. The orientation angle $\theta$ is defined as shown in Figure 7, and the averaged orientation angle is plotted as a function of the distance from the muscovite surface in Figure 8. Cosine of the angle was computed and averaged over time and molecules, which was then presented in terms of the angle $\theta$ (°), as calculated in previous studies [50,51,52]. If a water molecule rotates randomly, the time-averaged angle becomes 90°.

As shown in the snapshot in Figure 5a, hydrogen atoms of H₂O molecules adsorb directly on the muscovite, which results in the high orientation angle ($\theta$~150°) in the first layer (0.1 nm < d < 0.2 nm) in Figure 8.

In the second layer (0.2 nm < d < 0.5 nm), the profile of the orientation angle has two extreme values (~75° and ~115°). This indicates that the water molecules tend to exhibit a specific orientational order in the second layer, although not so as in the first layer.

In the third layer (0.5 nm < d < 0.7 nm), a smaller deviation from 90° is observed in Figure 8. This means that the third layer is less structured than the first and second layers.

The profiles at the ambient condition are in good agreement with the result of MD simulation by Van Lin et al. (2019) [22]. As the temperature increases, the deviation from 90° is gradually reduced (see an enlarged view in Figure S9 in Supplementary Materials). This indicates that the water molecules at higher temperature are more likely to rotate freely than those at lower temperature.

3.1.4. 2D Maps of the Orientation Angle

To gain more insights into the orientation angle $\theta$, the 2D maps of the orientation angle on the muscovite (0 0 1) surface are plotted in Figure 9. Additionally, in order to analyze the difference of the 2D profiles on K⁺-site and o-site, the 1D profiles of $\theta$ on K⁺-site and o-site are plotted in Figure 10. Moreover, the distributions of $cos\theta$ in each of the three layers are plotted in Figure 11.

• The first layer

High-angle areas are observed at o-site. This trend is also seen in the 2D density map of oxygens and hydrogens in the first layer, as shown in Figure 4. Figure 10 shows that high $\theta$ regions in the first layer come from the orientation on o-site. This is interpreted from the 1D density profiles in Figure 6; adsorption peaks are seen in the order of hydrogen and oxygen.

Additionally, the distribution of $cos\theta$ exhibits a sharp peak at $\theta$~180° in Figure 11a. By comparing the distributions at 298 K, 473 K, and 673 K at 50 MPa, we see that the peak height is reduced as the temperature increases. This shows that they can more easily rotate at high temperatures, though the interaction from the muscovite surface constrains H₂O molecules to some extent.

• The second layer

In the second layer, the distribution profiles on K⁺-site and o-site are clearly different. The former is higher than 90°, while the latter is lower than 90°.

These different patterns can be explained by the 1D profile of $\theta$ on K⁺-site and o-site in Figure 10. In the second layer, $\theta$ is higher than 90° while it is 40°–50° on K⁺-site. This leads to the different patterns in the 2D maps. Different orientations on K⁺-site and o-site can be understood by looking into the 1D density profile in Figure 6. On o-site, the adsorption peaks are seen in the order of hydrogen and oxygen from the muscovite surface. This indicates that H₂O molecules adsorb with high orientation angles. In contrast, on K⁺-site, the adsorption peaks are seen in the order of oxygen and hydrogen from the muscovite surface. This indicates that H₂O molecules adsorb with low orientation angles, directing the oxygen atoms toward the muscovite surface. In Figure 10b,c some outliers are observed at d~0.2 nm at high temperatures. This is probably because a small number of H₂O molecules can exist on K⁺-site due to the high kinetic energy of the H₂O molecules.

Temperature effect on the orientation angle in the second layer is manifest in the oscillatory profile for o-site in Figure 10 and its enlarged view in Figure S10 in Supplementary Materials. From these figures, the deviation from 90° around d = 0.35 nm is smaller as the temperature increases from 298 K to 673 K while that around d = 0.45 nm is larger in the order of 298 K > 473 K > 673 K. The large deviation is likely to be related to the peaks of the $cos\theta$ distribution at $cos\theta$ > ~0.7 and $cos\theta$ < ~−1, as shown in Figure 11b and Figure S11 in Supplementary Materials. Further studies are required to interpret the changes in the amplitude and the period of the fluctuation.

• The third layer

High- or low-angle areas are not clearly observed. This is consistent with the density profile in Figure 3 that shows that the third peak is less distinct than the first and second peaks. This indicates that water molecules are less structured in the third layer. As shown in Figure 11c, the distribution is flatter than those for the first and the second layers. This means that the distribution gets closer to that in bulk phase seen in Figure 11d, which is almost a flat distribution.

3.1.5. Hydrogen Bonding

Hydrogen bonding is considered to play a key role in forming the hydrogen network structure of water, not only in bulk water [53,54], ice [55], and hydrate [56], but also in adsorbed water [24]. Here, we focus on the number of hydrogen bonds per H₂O molecule to discuss the change in the structure of the adsorption layer due to the change in temperature. We calculated based on the existing criteria [57] (see Figure S12 in Supplementary Materials), which is commonly used with “gmx hbond” command in GROMACS [44]. In Figure 12, is plotted as a function of distance from muscovite surface. Figure 12a–c shows between H₂O-H₂O, H₂O-muscovite, and their sum, respectively. Compared to the result at 0.1 MPa/298 K with a previous MD simulation study by Wang et al. 2005 [24] in which is calculated for water/muscovite interfacial system, all of between H₂O-H₂O, H₂O-muscovite, and summation of them at the ambient condition are in good agreement with their result. Additionally, we compared of H₂O-H₂O in the bulk region for verification of the temperature effect. A number of studies have measured and calculated of H₂O-H₂O in a wide range of pressure and temperature conditions, and our results are in good agreement with them [58,59,60].

Where the distance from the surface is smaller than 0.5 nm, of H₂O-H₂O decreases as the distance gets smaller. This is consistent with previous studies using water/hydrophilic surfaces [24,61,62] as well as using hydrophobic surfaces [63,64]. Zhang et al. 2009 reported that the interactions between hydrophilic surface and water frustrate the hydrogen bonding network of water, resulting in the reduction of of H₂O-H₂O [61]. Considering the strong hydrophilicity of muscovite, the decrease of near the surface in Figure 12a is probably caused by the same mechanism. There is a small hump at around 0.2 nm from the surface. It gets smaller as the temperature increases, and of H₂O-muscovite increases as the distance gets smaller up to around 2. The temperature effect on of H₂O-muscovite is not visible as seen in Figure 12b. As a result, total is the highest at the first layer and it decreases at the second layer, with the lowest value around 0.3 nm from the surface. Our analysis suggests that the third layer is almost the same as the bulk region from the viewpoint of the hydrogen bonding as well as the orientation of H₂O molecules. These results clearly suggest that the destructuralization of the water adsorption layer is affected by the decrease in the number of hydrogen bonding.

In this study, we chose the SPC/E model as the force field for H₂O for consistency with the CLAYFF for muscovite in the same way as in previous studies [22,23,27,28,29,30,31]. It is, however, reported that the SPC/E model overestimates the energy of hydrogen bonding of H₂O molecules [65]. Therefore, further investigation into the effect of choice of the force field on the number of hydrogen bonding is necessary in a future work, especially for more precisely evaluating the contribution of hydrogen bonding to the water structure in the adsorption layers.

3.2. Dynamic Properties

One-dimensional self-diffusion coefficient $D_f$ is calculated from the mean square displacement as follows:

(2)

where $d_j\left(t_0\right)$ and $d_j\left(t\right)$ are the distance between a molecule $j$ and the surface at timestep $t_0$ and $t$. $N$ is the number of H₂O molecules, and < > indicates an average over $t_0$.

$D_f’$ for a molecule $j$ is calculated from:

(3)

where < > indicates an average over $t_0$. $D_f$ as a function of $d$ is then evaluated using $\rho \left(j,d\right)$, which is the probability to find molecule $j$ at $d$ [66]:

(4)$D_f\left(d\right)=\frac{{{\displaystyle \sum }}_{j=1}^N{D}_f’\left(j\right)\rho \left(j,d\right)}{{{\displaystyle \sum }}_{j=1}^N\rho \left(j,d\right)}$

$D^{*}{}_f\left(d\right)$, which is the calculated self-diffusion coefficient of water molecules normalized by the bulk value, is plotted in Figure 13. $D_f\left(d\right)$, without normalization, is plotted in Figure S13 in Supplementary Materials. $D_f\left(d\right)$ at bulk regions is in good agreement with the self-diffusion coefficient of water molecules obtained by the MD simulation in the previous study [67,68,69,70]. $D_f\left(d\right)$ at high temperature conditions (e.g., 673 K) is more than 10 times higher than that at ambient temperature. It is confirmed that diffusivity increases as the temperature increases, and the diffusivity tends to decrease as approaching the muscovite surface. The same trend was seen in a water/quartz interfacial system in a previous MD simulation [66]. The $D^{*}{}_f\left(d\right)$ profiles suggest that the interactions from the interface have some influence up to d~2.0 nm, which is far beyond the range of three layers, though the strength of influence depends on the pressure and temperature conditions. In all cases, $D^{*}{}_f\left(d\right)$ decreases as $d$ decreases. Focusing on the range of three layers, $D^{*}{}_f\left(d\right)$ at $d$ < 0.5 nm is less than 60% of that in the bulk region. $D^{*}{}_f\left(d\right)$, however, slightly increases in the very close vicinity of the surface because a repulsive force from the muscovite surface acts on the water molecules when they approach the vicinity of the muscovite surface.

At low temperature conditions such as 298 K and 323 K, $D^{*}{}_f\left(d\right)$ decreases sharply in the region closer to 0.5 nm from the muscovite surface. At 50 MPa/298 K, for example, the sharp peak of $D^{*}{}_f\left(d\right)$ in the first layer is 75% lower than that in the bulk region. Moreover, at these temperature conditions, a difference of $D^{*}{}_f\left(d\right)$ between d < 0.2 nm and 0.2 nm < d < 0.5 nm is clear; therefore, this result shows that dynamic properties are fairly different at each of the three layers. The first and second layers are much less diffusive than the third layer. At higher temperatures such as 423 K and 473 K, $D^{*}{}_f\left(d\right)$ decreases as d decreases, though the decrease in the first and the second layer is less clear than that in lower temperature conditions. At temperature conditions higher than 523 K, $D^{*}{}_f\left(d\right)$ in the first and the second layers is not clearly different.

4. Implication to Hydration Forces and Wettability

The layered ordering of the interfacial water is the origin of the hydration forces (also called solvation forces or structural forces) between hydrophilic surfaces in aqueous solution, which are essential to interpret a variety of phenomena such as colloidal dispersion or the swelling of clays [18,71]. For smooth surfaces such as the (0 0 1) plane of muscovite, the hydration forces consist of monotonic repulsive forces and oscillatory forces [71]. Uhlig et al. (2021) [23] confirmed that the oscillation is derived from the layered structure of water on muscovite by the AFM measurement of hydration forces of water/muscovite interface and MD simulation of the adsorption structure. Since hydration forces originate from the structure of the water adsorption layer, they are affected by temperature change that alters the structure of interfacial water as shown in our calculation. Structural forces derived from hydrogen bonding, which cannot be interpreted by the DLVO theory, largely contribute to hydration forces [18,72]; therefore, the decrease in the number of hydrogen bonding shown in 3.1.5 implies the decrease of the hydration forces. In fact, raising temperature is expected to decrease the structural repulsion [32,33], though this has not been measured or calculated at high-temperature conditions such as 673 K. Although our results imply that both magnitudes and periodicity of the fluctuation of the hydration forces are likely to change due to the temperature change, in order to understand quantitative effect it is necessary to interpret the changes in the density and the orientation profiles from the energetic viewpoint.

Hydration forces have a crucial role in determining wettability [72,73]. In previous studies, in water/gas/hydrophilic mineral or water/oil/hydrophilic mineral systems, adsorbed water film, the thickness of which is on the order of a nanometer, is formed between the non-wetting phase (gas or oil) and the mineral surface as shown in Figure 14. Mugele et al. (2015) [74] measured the thickness in decane/brine/muscovite systems, and it was from less than 1 nm to 10 nm depending on the ion species and concentration. MD simulation of decane/brine/muscovite systems shows that the thickness is from 0.5 nm to 1.3 nm depending on the salinity [75,76]. Based on the interfacial thermodynamics, the potential energy of the water film is known to be a crucial factor that determines the contact angle based on the augmented Young Laplace equation [72,73,77,78,79,80,81]. Churaev (1988) [32] pointed out the decrease in the structural repulsion as a cause of the worsening of the hydrophilicity reported by previous study [82]; however, quantitative analysis focusing on the changes in the water film due to the temperature change has not been conducted and it is expected as a future work.

The film thickness depends on a thermodynamic equilibrium condition, and there can be multiple thicknesses [72,73,80]. A previous study based on the augmented Young Laplace equation has shown that the two thicknesses $h$ (thinner: $h_1$= 0.3 nm, thicker: $h_2$ = 0.7 nm) have different contact angles for the CO₂/brine/muscovite system [80]. Hydration forces largely affect the stability of the film in this system, and these two thicknesses correspond to the second and third layers of the adsorbed water on the muscovite obtained in this study. Although the dynamic properties of these films were not known, from Figure 13, $D^{*}{}_f\left(d\right)$ of the third layer is 25%–40% lower than that of the bulk, and that of the second layer is 40%–75% lower than that of the bulk, which are highly temperature-dependent. Due to the difference in fluidity at the boundary between the second and third layers, the third layer will disappear easily depending on flow and thermodynamic conditions, while the second layer, that is, a thinner film, will remain on the surface. Although the second layer is more structured than the third layer, it is not perfectly immobile but has diffusivity. Therefore, transport of ions or some other molecules through the film can occur differently, depending on the film thickness and its property. This can be a key to understanding the time dependence of interfacial phenomena such as wettability alteration.

Future work is expected to study the property of the water adsorption layer under the existence of gas phases. We are going to work on the free energy calculation by following Morishita et al. [83,84] to study the stability of the adsorption layer. Moreover, calculation of the diffusivity at each adsorption layer in parallel to the mineral surface following the survival probability of molecules [85] at each layer is another future outlook for quantitative analysis of the mass transport in the thin water film.

5. Conclusions

We performed molecular dynamics (MD) simulations to study the structure and dynamics of water adsorption layers on muscovite mica at several pressures (0.1 MPa, 10 MPa, and 50 MPa) and temperatures (298 K, 323 K, 373 K, 423 K, 473 K, 523 K, 573 K, 623 K, and 673 K). We revealed structural changes of the adsorption layer based on analyses on the 1D and 2D profiles of density and orientation of H₂O. Additionally, we calculated self-diffusion coefficients of H₂O molecules in the adsorption region and the bulk region. The results show that the difference between these dynamical properties near the interface and in the bulk region is reduced as the temperature increases. Finally, we discussed implications to hydration forces of muscovite surface and wettability, which are significant interfacial properties in multiphase systems such as gas/brine/mineral systems, focusing on the structure and dynamics of the water adsorption film.

Supplementary Materials

The following are available online at https://www.mdpi.com/article/10.3390/w13091320/s1, Table S1: Distance between two muscovite surfaces, Figure S1: Normalized density profiles of H₂O, Figure S2: Normalized density profiles of oxygen of H₂O, Figure S3: Normalized density profiles of hydrogen of H₂O, Figure S4: Normalized density profiles of K⁺ on the muscovite surface, Figure S5: Normalized density profiles of K⁺, O (H₂O), and H (H₂O), Figure S6: Schematic images of the lattice of (0 0 1) surface of muscovite, Figure S7: 2D normalized density maps of K⁺ on the muscovite surface, Figure S8: Schematic image of one K⁺-site for calculating profiles of density and the orientation angle, Figure S9: Enlarged view of Figure 8c, Figure S10: Enlarged view of the part of Figure 10, Figure S11: Distributions of cosθ on K⁺-site and o-site at the second layer at 50 MPa, Figure S12: Schematic image of the definition of hydrogen bond, Figure S13: Self-diffusion coefficient of H₂O molecules vs distance from the muscovite surface.

Author Contributions

M.S. (Masashige Shiga) designed and performed the research; M.S. (Masashige Shiga), M.A., M.S. (Masao Sorai) and T.M. discussed the results and wrote the manuscript. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

The computation was carried out using the computer resource offered under the category of General Projects by Research Institute for Information Technology, Kyushu University. We thank anonymous reviewers for suggesting substantial improvements.

Conflicts of Interest

The authors declare no conflict of interest.

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

Figures and Table

Figure 1. A snapshot of the initial state of the calculation system (H2O/muscovite system). The area outside the blue line is the area duplicated by the periodic boundary conditions. Each atom is colored differently—oxygen: red; hydrogen: white; aluminum: cyan; silicon: yellow; potassium: orange.

Figure 2. Schematic image of muscovite structure (a) XY view (b,c) XZ views of the surface (d) K+-site (yellow: Si, blue: Al, red: O, white: H, orange: K).

Figure 3. Density profiles of H2O molecules vs. distance from the muscovite surface.

Figure 4. 2D normalized density maps of O and H atoms of H2O molecules at (a) the first layer, (b) the second layer, and (c) the third layer; the lattice on the map is from Figure S6 in Supplementary Materials.

Figure 5. Snapshots of the muscovite surface, K+, and H2O molecules.

Figure 6. Normalized density profiles of O, H, and K+ on K+-site and o-site.

Figure 7. Definition of the orientation angle of H2O.

Figure 8. Profiles of the orientation angle θ of H2O molecules from the muscovite surface at (a) 0.1 MPa (b) 10 MPa (c) 50 MPa; a snapshot of a H2O molecule at each of three extreme values of θ is shown.

Figure 9. 2D Maps of the orientation angle θ at 50 MPa; the lattice on the map is from Figure S6 in Supplementary Materials.

Figure 9. 2D Maps of the orientation angle θ at 50 MPa; the lattice on the map is from Figure S6 in Supplementary Materials.

Figure 10. Profiles of the orientation angle θ on K+-site (red) and o-site (green).

Figure 10. Profiles of the orientation angle θ on K+-site (red) and o-site (green).

Figure 11. Distributions of cosθ at (a) the first layer, (b) the second layer, (c) the third layer, and (d) bulk region (2.0 nm < d < 3.0 nm). Each pressure condition is 50 MPa.

Figure 12. The number of hydrogen bonding per H2O molecule vs distance from muscovite surface.

Figure 13. The normalized self-diffusion coefficient of H2O molecules D*fd vs distance from the muscovite surface.

Figure 14. Schematic image of the thin film with two different thicknesses (thinner: h1; thicker: h2) in the gas/water/mineral system.

Table 1

Pressure and temperature conditions.

Pressure (MPa)	Temperature (K)
0.1	298, 323, 373, 423, 473, 523, 573
10	298, 323, 373, 423, 473, 523, 573, 623
50	298, 323, 373, 423, 473, 523, 573, 623, 673