Machine-learning enhanced simulations predict graphene is hydrophobic and microscopically not wetting transparent

Domain: biology
Published: 2026-04-02
Editorial outcome: minor_revision

Abstract
The interaction between graphene and water is fundamental to applications ranging from filtration to nano-electronics, yet the intrinsic wettability of graphene remains a subject of longstanding debate. In particular, it is disputed whether graphene is wetting transparent—transmitting the wettability of an underlying substrate to the surface. Here we show using machine-learning enhanced molecular dynamics simulations that pristine graphene is intrinsically hydrophobic and microscopically not wetting transparent. By simulating vibrational sum-frequency generation spectra, we demonstrate that the apparent hydrophilic signatures often observed for monolayer graphene on hydrophilic substrates originate not from wetting transparency, but from signal cancellation caused by water molecules intercalated between the graphene and the substrate. Furthermore, we find that while water intercalation is thermodynamically favorable for monolayers, it becomes unfavorable for multilayer graphene, explaining experimentally observed thickness-dependent wetting behaviors. These findings provide a unified microscopic framework for understanding graphene-water interactions and clarify the crucial role of confined water in two-dimensional material interfaces.

Introduction
Since the successful mechanical exfoliation of graphene in 2004,[1] two-dimensional (2D) materials have garnered widespread scientific interest owing to their unique physical and chemical properties, which are fundamentally distinct from those of their bulk counterparts.[2–4] Graphene, a single layer of carbon atoms arranged in a hexagonal lattice with sp2 hybridization, is one of the most representative 2D materials.[5] It possesses a combination of remarkable properties, including high electrical and thermal conductivity, excellent optical transparency, superior chemical stability, outstanding mechanical strength, and complete impermeability to atoms and molecules.[6–11] These remarkable properties have enabled a wide range of applications, such as desalination of water, transparent conductive electrodes, electrical energy storage, and electrocatalysis.[12–16] Notably, many of these applications involve the interaction of graphene with water or aqueous electrolytes.[17–19] This interaction enforces a distinct molecular structure of water, expressed through orientation of water molecules and the hydrogen-bond network, both of which are crucial for the optimal performance of devices like proton exchange membrane fuel cells[20–22] and nanofluidic devices[23, 24].

One of the most widely used experimental techniques for evaluating graphene-water interactions is the measurement of the water contact angle (WCA).[25–27] However, WCA values are influenced by a variety of factors, including the presence of an underlying substrate,[28–30] the defect density on both the graphene and the substrate,[31–34] the number of graphene layers,[30,35] and surface contaminants[36–38]. These factors can result in a wide range of reported WCA values, ranging from 30◦ to 127◦.[39–41] Given that materials with WCA below 90◦ are considered hydrophilic and those above 90◦ hydrophobic, there remains ongoing debate as to whether graphene is intrinsically hydrophilic or hydrophobic.[42] Although WCA measurement is a useful tool for characterizing surface wettability, it is a macroscopic observable and does not directly reveal molecular-level details of the interfacial water structure. Therefore, a comprehensive understanding of the wettability of graphene requires molecular-scale insight into the graphene-water interactions that govern its observed WCA behavior.

Vibrational sum-frequency generation (vSFG) spectroscopy[43–47] is a powerful technique for probing the molecular structure of water at interfaces.[48–54] As a second-order non-linear optical process, vSFG requires breaking the inversion symmetry to generate a finite signal. This requirement makes the technique inherently surface-specific, as bulk water is centrosymmetric and thus does not contribute to the optical response. Consequently, the vibrational signals detected by vSFG arise exclusively from interfacial water molecules. Previous experimental and simulation results from Wang et al.[55], together with independent experimental measurements by Xu et al.[56], consistently demonstrated that pristine graphene floating on water exhibits hydrophobic behavior. More recently, Rashmi et al.[57] employed path-integral molecular dynamics with the MB-pol potential to further confirm the hydrophobic nature of neutral graphene suspended on water. In contrast, when supported on a calcium fluoride (CaF2) substrate, the measured vSFG spectra indicated a hydrophilic character of the graphene layer.[35, 54, 58] Given the hydrophilicity of the substrate and the atomic thickness of graphene, some researchers proposed that the wettability of the underlying substrate significantly influences the overall interfacial properties of graphene and water. This phenomenon was referred to as wetting transparency and was proposed to explain the observed wettability of the graphene surface.[30] However, the extent to which graphene is transparent to wetting remains a subject of ongoing debate. Graphene has been reported to be fully wetting-transparent,[30, 35, 59, 60] partially wetting-transparent,[37, 61–63] or fully wetting-opaque[41, 64–67].

In addition to substrate effects, intercalated water molecules may also influence the shape of vSFG spectra. Ohto et al.[68] simulated the vSFG spectra of pristine graphene and found that when all water molecules were located on one side of the graphene sheet, the spectrum closely resembled that of the water/air interface, consistent with the hydrophobic nature of graphene. However, the presence of even a small number of water molecules on the opposite side of the graphene layer significantly reduced the peak amplitude of the dangling O-H bond, indicating a diminished hydrophobic response. In experimental work, Temmen et al.[69] used non-contact atomic force microscopy and Kelvin probe force microscopy to characterize the thickness of graphene sheets and the intercalated water layer. They found that even after annealing to 750 K, it was nearly impossible to completely remove the water molecules trapped between graphene and the substrate CaF2. Li et al.[70] investigated the behavior of intercalated water by calculating the migration energy barrier of a single water molecule confined between monolayer graphene and substrate CaF2, using the climbing image nudged elastic band method. They reported a migration barrier of approximately 0.19 eV. They also calculated the diffusion coefficient of water at 300 K, finding it to be 8.35 × 10−6 cm2/s. These results, a low migration barrier, a high diffusion rate, and the persistence of water even after high-temperature annealing, collectively suggest the presence of intercalated water between graphene and the CaF2 substrate.

Beyond substrate effects and intercalated water, the number of graphene layers can also influence the experimentally measured vSFG spectra. Kim et al.[35] systematically measured the WCA and vSFG spectra of multilayer graphene. Their vSFG measurements revealed that when the number of graphene layers exceeded four, the vSFG spectra exhibited pronounced changes. These changes were attributed to the transition of the graphene film from wetting transparent to translucent and finally to opaque as the number of layers increases from one to six. However, the fundamental mechanisms underlying these dramatic changes in the vSFG response with increasing graphene thickness were not fully understood.

Although experimental vSFG provides valuable information about interfacial water structure, it remains challenging to control atomically the interface to disentangle the various factors, such as graphene-water interactions, graphene layer effects, substrate influence, and water penetration. To achieve a comprehensive microscopic understanding, atomic-scale simulations are indispensable. Ab Initio molecular dynamics (AIMD) simulations can provide such insights with quantum-chemical accuracy that allows for a reliable description of the water orientation and hydrogen-bond network at the graphene/water interface. Unfortunately, AIMD simulations are generally limited by high computational cost and insufficient sampling. Accurate simulation of vSFG spectra often requires nanosecond time-scale trajectories to ensure statistical convergence.[71, 72] To solve this issue, machine learning interatomic potentials (MLIP) have recently become impressively successful in interpolating the quantum-chemical potential energy surface, reaching AIMD-level accuracy while dramatically improving computational efficiency.

Accurate simulation of the water structure at the graphene/water interface is a complex task. The orientation-dependent adsorption energy of single water on graphene varies by less than 0.043 eV (1 kcal/mol), necessitating simulations that exceed the threshold of chemical accuracy.[73–77] Moreover, the graphene-water interactions are modulated by several factors, including the substrate,[30, 78] the number of graphene layers,[35, 79] the defect of graphene,[34] the thickness of the water slab,[80, 81], external electrical field,[56, 57, 82, 83] and the presence of intercalated water molecules between the graphene and the substrate.[84–86] Capturing such a complex, multi-element, and many-atom system with sub-chemical accuracy pushes the capabilities of MLIP beyond what has been explored in the existing literature. In addition, previous theoretical studies of vSFG spectra at the graphene/water interface have primarily focused on simplified systems, such as monolayer graphene with a few water molecules,[55, 68, 87, 88] or systems with the CaF2 substrate and a monolayer graphene containing a small number of intercalated water molecules.[86].

In contrast to earlier experimental and theoretical studies—typically restricted to either suspended graphene/water interfaces or supported monolayer graphene under fixed conditions, we developed chemically accurate MLIPs based on the Atomic Cluster Expansion (ACE) method,[89, 90] trained in extensive density functional theory (DFT) datasets. These MLIPs enable long-time, large-scale molecular dynamics (MD) simulations of graphene/water interfaces spanning a broad range of conditions, including different numbers of graphene layers, the presence or absence of a CaF2 substrate, and varying amounts of intercalated water. The resulting MD trajectories allow us to systematically investigate the orientation of water molecules at the graphene/water interface and to simulate the corresponding vSFG spectra.

In this work, we show the individual roles of the substrate, graphene layer number, and intercalated water in governing the vSFG response by employing machine-learning enhanced MD simulations. We identify signal cancellation between intercalated and adsorbed water as the microscopic origin of the hydrophilic vSFG features observed for supported monolayer graphene on the CaF2 substrate, challenging the hypothesis of wetting transparency. Furthermore, we quantify how the thermodynamic driving force for water intercalation varies with graphene thickness, explaining the experimentally reported thickness-dependent transition in wettability and confirming that pristine graphene remains intrinsically hydrophobic regardless of the layer number.

Results
2.1 Atomistic models

To systematically capture substrate, water intercalation, and graphene layer effects (Fig.1a), we constructed a series of atomistic models and adopted a consistent naming scheme. Each system is labeled as SnLm, where n denotes the number of intercalated water molecules in our simulation cell and m indicates the number of graphene layers. The water slab, which is present in all models, is not explicitly included in the naming scheme. If S, n, or L is omitted in the name, it implies the absence of the substrate, intercalated water, or graphene layers, respectively. According to this convention, L1 and L5 denote systems comprising one- and five-layer graphene, respectively, without substrate and intercalated water, but with a water slab placed above the graphene. The label S refers to the bare CaF2 substrate with a water slab. For supported systems without intercalated water, SL1, SL4, and SL5 represent the substrate combined with one-, four-, and five-layer graphene, respectively, along with an overlying water slab. Finally, SnL1 and SnL4 denote systems with one- and four-layer graphene, respectively, that contain n intercalated water molecules confined between the graphene and the substrate, in addition to the overlying water slab.

Li et al.[70] compared the binding energies of various phases of monolayer ice and monolayer amorphous water confined between monolayer graphene and the CaF2 substrate. Their results showed that one particular ice phase was more stable than the other configurations. Motivated by their findings, we extended our study beyond models with discrete numbers of intercalated water molecules to include several ordered monolayer ice structures, as well as one double-layer amorphous water confined between the substrate and monolayer graphene. The specific structural configurations are illustrated in Suppl. Fig. 1. According to our naming convention, SHL1 denotes a monolayer of hexagonal ice (ice-Ih), while SPL1 corresponds to a monolayer of pentagonal ice. Furthermore, SS1L1 and SS2L1 represent two different monolayer square ice structures, and SDL1 refers to a double-layer amorphous water configuration. This naming scheme enables a systematic comparison across different structural and interfacial configurations.

2.2 Potential training and accuracy

To train the MLIP, a representative dataset was first constructed using on-the-fly machine learning force field MD simulations.[91–93] This approach significantly accelerates MD simulations for complex systems such as SL5 and SnL4. Atomic configurations, along with their corresponding DFT computed energies and forces, were extracted from the MD trajectories to form the initial dataset. The ACE method[89, 90] was then employed to train the MLIPs. Given the structural complexity and diversity of systems considered (see section Atomistic models), it is impractical to develop a single potential capable of accurately describing all configurations, since the ACE model is inherently local and short-ranged. Merging all training datasets to build a single universal potential would require a dramatically larger basis and/or body order to reach the same level of precision. In particular, the relevant water–water separations differ strongly between SnL1 and SnL4 systems: intercalated and slab waters fall inside the cutoff in SnL1, but are separated by ≳15 Å in SnL4, i.e., far beyond the cutoff range of the short ranged expansion. A universal short-range potential would therefore have to interpolate between physically distinct regimes, which we expect to degrade its accuracy for both.

Considering that our simulations demand highly accurate interatomic potentials, system-specific ACE models provide the most reliable approach. Accordingly, multiple system-specific ACE potentials were trained: ACE-L1 for L1, ACE-L5 for L5, ACE-S for S, ACE-SL1 for SL1, ACE-SL4 for SL4, ACE-SL5 for SL5, ACE-SnL1 for all monolayer graphene intercalated water systems, and ACE-SnL4 for all four-layer graphene intercalated water systems. The system specific potentials are not intended to be transferable across distinct systems. For example, applying ACE-L1 to the SnL4 system causes the MD simulation to terminate within a few steps, as the extrapolation grade value becomes prohibitively large, indicating that the configuration lies far outside the training domain and that the potential cannot guarantee reliable predictions.

To ensure that each potential is reliable for large-scale and long-timescale MD simulations, we further refined MLIPs using an active learning (AL) strategy.[94] Fig. 1c schematically illustrates the methodological workflow, from on-the-fly machine learning force field MD simulations to iterative MLIP training and refinement, followed by large-scale, long-timescale MD simulations and subsequent statistical and spectroscopic analyses. Detailed descriptions of each step are provided in the Methods section. The accuracy of the trained potentials was evaluated using the root mean square error (RMSE) of energies and forces on the test dataset, as summarized in Suppl. Table 1. All trained MLIPs exhibit energy RMSE values no more than 0.50 meV/atom and force RMSE values below 40 meV/Å, indicating high accuracy of our MLIPs.

To ensure the accuracy and stability of the large-scale, long-timescale MD simulations, we continuously monitored the extrapolation grade value during both the equilibration and production MD simulation stages. The extrapolation grade values for the production MD simulation are shown in Suppl. Fig. 11-29. The time evolution of the temperature of all water molecules, along with the total internal energy of each system during the NVE ensemble MD simulations, is presented in Suppl. Figs. 30-48. Across all NVE ensemble MD simulations, the total energy of each system fluctuates by approximately 0.1 eV, further confirming the stability and reliability of the trained MLIPs. The atomic probability distributions, water density profiles along the z-axis and radial distribution functions of gOO(r) and gOH(r) for all systems are presented in Suppl. Figs. 49 - 67, Suppl. Figs. 68 - 86, and Suppl. Figs. 87 - 90, respectively.

2.3 Effect of various graphene layers

It is widely accepted that graphite exhibits hydrophobic properties. However, the wettability of monolayer and multilayer graphene remains a topic of debate.[42] In this study, we first considered monolayer (L1) and five-layer graphene (L5) systems without substrates and defects. To analyze the orientation of interfacial water molecules, we defined two angles (Fig. 1b): σ, which describes the orientation of the lower coordination O-H bond with respect to the Z-axis, and ϕ, representing the angle between the normal vector of the water plane and the Z-axis. As shown in Fig. 2b and 2e, for the first water layer adjacent to graphene, both L1 and L5 systems exhibit three distinct types of water molecule orientations. Among them, the water molecules closest to the graphene surface often have one O-H bond pointing directly toward the surface (also known as dangling O-H). Notably, the L5 system shows a significantly higher number of such dangling O-H bonds compared to the L1 system. This suggests that L5 is more hydrophobic than L1, consistent with experimental measurements of contact angles, which showed an increase in hydrophobicity with the number of graphene layers.[30, 35, 79] In contrast, the angular distribution at the water/air interface does not exhibit such well-defined orientation types (Suppl. Fig. 91), reflecting the weaker constraints imposed by the gas/liquid interface compared with the stronger confinement at the solid/liquid interface.[87]

We then calculated the vSFG spectra for L1 and L5 systems, using the well-studied water/air interface as a reference. Experimentally, the vSFG peak corresponding to the dangling O–H stretch of water/air interface appears at approximately 3700 cm−1.[56, 95–98] As shown in Fig. 2c, the simulated imaginary part, Im(χ(2)), of the vSFG spectrum of water/air interface shows the dangling O–H peak at 3730 cm−1. The slight blue shift relative to experimental value is attributed to the combined effect of classical nuclei (neglecting nuclear quantum effects) and the underlying density functional, which together are known to shift OH-stretch frequencies while largely preserving the overall spectral shape.[99, 100] The broad positive shoulder around 3600 cm−1 observed at the water/air interface has been attributed in previous studies to intramolecular vibrational coupling between O-H chromophores.[101] This feature is suppressed in the uncoupled case of dilute HOD in H2O. For L1, our simulated vSFG spectrum agrees well with previous simulations and experimental measurements.[55, 56, 87, 88]. Such consistency is further validated by the recent study of Rashmi et al.[57], employing high-accuracy path-integral molecular dynamics, which independently confirms the hydrophobic nature of the suspended monolayer graphene. The dangling O-H peak of L1 is red-shifted by 65 cm−1 relative to the water/air interface, which is attributed to O-H-π interactions that weaken the O-H bond and lower its vibrational frequency.[35, 56, 84] For L5, water molecules close to the graphene surface have a higher ratio of dangling O-H bonds and consequently less hydrogen-bonded O-H. This molecular arrangement is reflected in Im(χ(2)) spectra. Compared to Im(χ(2)) spectrum of L1, Im(χ(2)) spectrum of L5 shows a reduced amplitude around 3200 cm−1 (associated with hydrogen-bonded O-H stretch vibrations) and an enhanced amplitude near 3665 cm−1 (characteristic of the dangling O-H groups). These differences are more obvious in the square modulus |χ(2)|2 of the vSFG spectra (sum of the squared imaginary and real parts)(Fig. 2f). These findings further confirm that both monolayer and multilayer pristine graphene are intrinsically hydrophobic, with the degree of hydrophobicity increasing with the number of layers. This trend is consistent with the enhanced exposure of the dangling O-H groups and is supported by experimental contact angle measurements.[30, 35]

2.4 Substrate effect

Most experimental studies reporting vSFG spectra of graphene used CaF2 as a substrate, owing to its wide optical transparency range and the chemically inert nature of the fluorine-terminated (111) surface, which is stable in air, has low solubility in water, and lacks dangling bonds or charge traps.[35, 49, 58, 82, 84–86, 102–104] The vSFG spectrum of monolayer graphene on the CaF2 substrate often exhibits features indicative of hydrophilic behavior.[35, 54, 58] Since CaF2 is well known to be hydrophilic and has a strong affinity for water molecules,[105–107] this observation appears to contradict the intrinsic hydrophobicity of pristine monolayer graphene, as confirmed by both simulations and experimental measurements.[55–57, 68] To explain this apparent contradiction, some investigators have proposed the concept of wetting transparency,[30] where the wetting properties of the underlying substrate (CaF2) dominate the observed interfacial behaviour despite the presence of the graphene layer.

To investigate the substrate effect, we computed water orientation with respect to the substrate or graphene surface, as well as vSFG spectra, for several distinct systems. These include the bare CaF2 substrate with a 512 water molecule slab (S), shown in Fig. 3a, and the monolayer graphene supported by CaF2 substrate with a water slab (SL1) as depicted in Fig. 3d. We also examined multilayer systems, specifically the CaF2 substrate with four-layer graphene and the water slab (SL4), and the CaF2 substrate with five-layer graphene and the water slab (SL5), shown in Suppl. Fig. 94a and Suppl. Fig. 94d, respectively.

Based on the atomic distribution analysis of system S (Suppl. Fig. 51), the average position of the first water layer is closer to CaF2 surface than to the graphene surface (SL1, Suppl. Fig. 52). This proximity is attributed to the stronger interactions between water oxygen atoms and calcium ions on the CaF2 surface, in contrast to the relatively weak physisorption observed between water molecules and graphene. As a result, the fraction of dangling O-H bonds near the CaF2 surface is significantly reduced, reaching only about 60% of that observed in the L1 or L5 systems. The exact number of dangling O-H bonds is provided in Suppl. Note 10. In contrast to graphene-based systems (L1 and L5), the second type of interfacial water in S (indicated by the red arrow) is oriented toward the bulk water rather than lying parallel to the surface. This distinct interfacial arrangement leads to a markedly different Im(χ(2)) spectral profile (Fig. 3c). Unlike graphene systems (L1 and L5), the CaF2 and water interface does not exhibit a positive peak in the Im(χ(2)) spectrum around 3665 cm−1, which in vSFG spectra is indicative of dangling O–H groups. Given the sensitivity of vSFG to molecular orientation, the presence of even a small number of dangling O-H bonds would generally produce this positive peak. The absence of such a feature in system S can be attributed to two main factors. First, there is a significantly reduced number of dangling O-H bonds. Second, the dangling O–H signal is partially cancelled by the non-hydrogen-bonded, second-type water molecules oriented toward the bulk. Moreover, our simulated vSFG spectrum for system S shows good agreement with previous simulation and experimental measurements,[35, 58, 107] showing only a single broad negative peak characteristic of hydrogen-bonded O-H. This consistency supports the validity of our structural model and interpretation.

The angular distribution of water molecules in SL1 (Fig. 3e), SL4 (Suppl. Fig. 94b), and SL5 (Suppl. Fig. 94e) systems is quite similar, and closely resembles that of the corresponding systems without a substrate (L1 and L5). In the first water layer adjacent to the graphene surface, three distinct configurations of water molecules are observed. Among these, the water molecules closest to the graphene typically adopt a configuration in which one O–H bond points toward the graphene surface. In the substrate-supported systems (SL1, SL4, and SL5), the number of dangling O-H bonds increases with the number of graphene layers and is consistently higher than the corresponding substrate-free systems (L1 and L5). The elevated number of dangling O–H bonds in the first water layer of substrate-supported systems suggests that the presence of the CaF2 substrate may enhance the hydrophobic character of the interface. This observation contrasts with the concept of wetting transparency.

In addition to the similar angular distributions of water molecules observed in systems L1, L5, SL1, SL4, and SL5, Im(χ(2)) spectra exhibit a generally consistent spectral profile across all these five systems. Although there are minor differences in the intensities of the positive and negative peaks, attributable to differing amounts of hydrogen-bonded versus dangling O-H groups, the overall shape of the spectra remains largely unchanged. Furthermore, the square modulus of the vSFG spectra, |χ(2)|2 (Fig. 3f), derived from our simulations, differs from the experimental spectra reported by Kim et al. [35]. Notably, no significant change is observed at the frequency of 3665 cm−1 as the number of graphene layers increases. These simulation results provide no evidence supporting wetting transparency in monolayer and multilayer graphene systems, suggesting that the substrate alone is insufficient to induce a substantial change in the vSFG spectral response.

2.5 Intercalated water effect

To rationalize the deviations between our simulated and experimental vSFG spectra for supported multilayer graphene, we first examine the possible role of water intercalation. Kim et al.[35] reported that the defect density (defects per carbon atom) in graphene with values of (13.6 ± 3.8) × 10−6 for monolayer graphene, (9.3 ± 2.6) × 10−6 for four-layer graphene, and (5.1 ± 1.5) × 10−6 for five-layer graphene. Given these extremely low defect densities, the graphene sheets in our simulation models were assumed to be defect-free. It is well established that pristine graphene is impermeable to water molecules.[108] Consequently, for multilayer pristine graphene, water cannot intercalate between the graphene layers from the water-facing side. Arif et al.[109] used atomic force microscopy to measure the thickness of multilayer graphene under relative humidities ranging from 10% to 80% and observed only minimal changes in thickness, concluding that water intercalation into multilayer graphene is negligible. Because our simulation cells employ periodic boundary conditions and thus contain no explicit edges, they represent the interior, defect-free, edge-sealed domains of multilayer graphene. Accordingly, and consistent with experimental observations, water intercalation between graphene layers on the water-facing side was not considered in our simulations. Instead, we only account for possible water intercalation between the graphene and CaF2 substrate.

Experimentally, no dangling O-H peaks are observed in the vSFG spectra for graphene films with up to three layers.[35] In contrast, our simulations for the monolayer graphene system exhibit a pronounced dangling O-H peak (Fig. 3f), suggesting that additional factors, such as intercalated water, may substantially influence the spectral features. It is well established that a thin layer of water readily forms on hydrophilic surfaces under ambient conditions.[110–113] When graphene is transferred onto a hydrophilic substrate, water molecules can be trapped beneath the impermeable graphene layer.[69, 114–116] Furthermore, first principles calculations for the CaF2(111)/graphene system have shown that once a water molecule is confined between monolayer graphene and the CaF2(111) substrate, its lateral diffusion between neighbouring adsorption sites proceeds with a very low migration barrier of approximately 0.19 eV, leading to a diffusion coefficient of about 8.35 × 10−6 cm2/s at 300 K.[70] Such facile diffusion implies that, once intercalated, water molecules can rapidly spread under the graphene layer, forming an extended confined layer. Additionally, during vSFG spectral measurements, where samples are typically immersed in water,[35, 58, 82, 84] capillary effect can further promote the intercalation of water molecules between the graphene and hydrophilic substrate. In the experimental work by Lee et al.,[117] who immersed annealed SiO2 supported graphene in water, real-time Raman spectroscopy was used to visualize the diffusion of water underneath the graphene layer. They found that the diffusion rate of water depends strongly on the hydrophilicity of the substrate. Motivated by these findings, we explicitly consider the effect of varying amounts of intercalated water molecules on the vSFG spectra.

In our model, a full monolayer of intercalated ice-Ih between the CaF2 substrate and monolayer graphene consists of 96 water molecules. In addition to this full monolayer, we also considered systems with partial coverage corresponding to one-quarter, one-half, two-thirds, and five-sixths of a monolayer of ice-Ih, containing 24, 48, 64, and 80 water molecules, respectively. These systems are denoted as S24L1 (Fig. 4a), S48L1 (Suppl. Fig. 95d), S64L1 (Fig. 4d) and S80L1 (Suppl. Fig. 96d), respectively. The system with a complete ice-Ih monolayer is labeled SHL1, as the water molecules form hexagonal rings (hexamers)(Suppl. Fig. 1a) stabilized by hydrogen bonding, characteristic of the ice-Ih structure. Beyond the monolayer ice-Ih configuration, we also explored alternative initial interfacial water structures, including different phases of monolayer ice and a bilayer of amorphous water, to include the possibility of metastable water structures that could influence wettability of graphene. The configuration of these water structures is shown in Suppl. Fig. 1.

We first considered systems containing intercalated water molecules in quantities equal to or less than that of a complete monolayer of ice-Ih. We generally found most of the intercalated water lying flat on the substrate for small water coverage (S24L1, Fig. 4a). When water coverage approached a full monolayer, water was increasingly oriented with O-H towards the substrate (Fig. 4e), becoming the dominating configuration at S80L1 (Suppl. Fig. 96f). For a complete monolayer coverage (SHL1)(Suppl. Fig. 97a), additionally, one more water configuration with O-H pointing towards the graphene surface is found.

It is noteworthy that the ideal hexagonal ring pattern characteristic of crystalline ice-Ih is not preserved after equilibration in our simulations (Suppl. Fig. 102g). This observation aligns with experimental observations from room-temperature atomic force microscopy, which also report the absence of well-ordered ice-Ih under similar conditions.[118] Specifically, while ice-Ih has been detected between graphene and substrates BaF2(111), it is not observed between graphene and substrate CaF2(111). This difference is primarily attributed to the lattice mismatch between the CaF2(111) surface and the crystalline structure of ice-Ih. These results suggest that the orientation and structure of intercalated water molecules are mainly governed by the degree of water coverage rather than a specific structural arrangement. It is also important to note that at sub-monolayer coverages, water does not form a homogeneous structure but instead organizes into discrete, localized clusters (Suppl. Fig. 102d, 102e and 102f).

We further examined how intercalated water modulates the vSFG spectral response. As shown in Fig. 4c and Suppl. Fig. 102a, the positive feature near 3665 cm−1 in the Im(χ(2)) spectra diminishes progressively with increasing intercalated water content. In the SL1 system lacking intercalated water, a shoulder at around 3500 cm−1 (Fig. 3c) arises from the antisymmetric O–H stretch of interfacial water molecules.[101] As additional water intercalate, the shoulder evolves into a distinct peak that shifts to more negative values. When isolating the contribution from the intercalated water molecules (Suppl. Fig. 102c), the Im(χ(2)) spectra reveal a distinct negative peak around 3650 cm−1 accompanied by a positive one near 3285 cm−1. The increasing amplitude of the negative feature reflects the growing population of water molecules with one O–H bond oriented toward the graphene, consistent with the dominant contribution observed in the SHL1 system. In experiments by Wang et al.[84, 85], who probed the effects of the substrate and intercalated water under varying external potentials. At high negative potentials, a negative peak near 3630 cm−1 was observed and attributed to the O–H stretching mode of CaOH species, whereas this feature disappeared at weakly negative or positive potentials. Although our simulations were conducted without any applied potential, increasing intercalated water similarly converts the dangling O–H positive feature into a negative peak, indicating that water intercalation alone is sufficient to induce spectral inversion. Furthermore, in SHL1, contributions from water molecules on both sides of graphene lead to a more complex Im(χ(2)) spectrum with three negative peaks(Fig. 4c). In the square modulus |χ(2)|2 of the vSFG spectra (Fig. 4f), the intensity of the peak at 3650 cm−1 first decreases and then increases with the addition of intercalated water. Collectively, these results suggest that the apparent hydrophilic signatures observed in the vSFG spectra of monolayer graphene on CaF2, despite the intrinsic hydrophobicity of pristine graphene, are likely due to the presence of intercalated water. During vSFG measurements, in which samples are immersed in water, capillary forces can drive water beneath the graphene, restructuring the interfacial water and modifying the vibrational response.

To gain deeper insight into how intercalated water structure and coverage affect the vSFG response, we further explored a range of crystalline monolayer ice models along with a representative double-layer amorphous water configuration. Across all monolayer ice systems, the initial crystalline order was lost during equilibration (Suppl. Fig. 103). Importantly, similar water orientations and vSFG spectral characteristics were observed in systems with identical surface coverage (Suppl. Fig. 103), independent of their initial structure. In contrast, the double-layer amorphous water configuration exhibited distinct spectral signatures, characterized by vSFG features. These results indicate that the vSFG response of intercalated water is more strongly influenced by intercalated water surface coverage than by the initial ice phase. Detailed structural and spectral analyses are provided in the Suppl. Note 11.

In addition to investigating monolayer graphene with intercalated water molecules, we further examined the influence of multilayer graphene on the behavior of intercalated water molecules and the corresponding vSFG spectra. Specifically, we selected a four-layer graphene system as the representative case to analyze the angular distribution and vSFG spectral features. The simulated systems consist of a CaF2 substrate, a four-layer graphene sheet, a water slab containing 512 water molecules, and varying amounts of intercalated water (24, 48, 64, and 96 molecules). These systems are referred to as S24L4 (Suppl. Fig. 100a), S48L4 (Suppl. Fig. 100d), S64L4 (Fig. 5a), and SHL4 (Suppl. Fig. 101d), respectively.

The angular distribution of intercalated water molecules between the CaF2 substrate and four-layer graphene exhibited a trend similar to that observed in systems with monolayer graphene. Specifically, as the amount of intercalated water increases, the orientation of water molecules evolves in a comparable manner to that of monolayer graphene. In Im(χ(2)) (Fig. 5b), taking S64L4 as an example, compared to the monolayer graphene systems(S64L1), the broad peak at the frequency of 2800-3500 cm−1 exhibits a blue shift. Additionally, the sharp peak near 3665 shows a higher amplitude in the four-layer graphene system. For thicker graphene layers, the reduced number of hydrogen-bonded O–H groups and the increased presence of dangling O–H bonds at the interface jointly lead to an interfacial arrangement consistent with a more hydrophobic character. The square modulus vSFG spectra, |χ(2)|2 (Fig. 5c), also follow a trend similar to that of the monolayer systems, but with a blue-shifted broad peak and an unobvious peak around 3650 cm−1. These spectral features further confirm the stronger hydrophobic character of multilayer graphene and its impact on the interfacial water structure.

2.6 Driving force for water intercalation

Although introducing an equivalent amount of intercalated water molecules results in similar spectral features in the simulated vSFG spectra of monolayer and four-layer graphene systems near 3665 cm−1 (Fig. 5c), experimental vSFG measurements exhibit pronounced differences between these two systems in this frequency range.[35] Independently measured vSFG spectra from different research groups for monolayer graphene on a CaF2 substrate do not exhibit a distinct peak near 3665 cm−1.[35, 54, 58] Although a weak shoulder feature around 3650 cm−1 was observed in the vSFG spectra measured by Dreier et al.,[58] this suggests that the monolayer graphene supported on CaF2 exhibits a hydrophilic character. In contrast, for graphene films thicker than three layers, a distinct peak (or shoulder) emerges in the frequency range of 3600-3700 cm−1,[35] whereas this feature is absent in systems with one to three graphene layers.

To understand this discrepancy, we calculated the free energy of intercalated water molecules and water molecules in the slab along the trajectory from internal energy data and an additional entropic contribution ( see Section 4.6 for details). We used the free energy per water molecule in each slab as the reference, setting its value to zero (Suppl. Table 8). Our results show that the free energy of intercalated water molecules is negative in the monolayer graphene system (Fig. 6a), indicating that water intercalation between graphene and the CaF2 substrate is thermodynamically favorable. Furthermore, the free energy remains nearly unchanged as the number of intercalated water molecules increases up to 12 per nm2 (S64L1 and S80L1). Beyond this coverage, further intercalation becomes energetically less favorable, as observed in systems such as SS1L1 (12.5 water/nm2), SS2L1 (12.5 water/nm2), and especially SDL1 (23 water/nm2). In contrast, the free energy of an intercalated water molecule becomes positive (> 1.5 eV) for the four-layer graphene system (Fig. 6a), indicating an energetically unfavorable intercalation process. At room temperature, a positive free energy penalty of about 1.7-2.1 eV per intercalated water molecule (Suppl. Table 8) corresponds to an exponentially small equilibrium population of confined water under four-layer graphene, even though the local lateral diffusion barrier within the CaF2(111)/graphene gap is expected to remain low. Thus, for multilayer graphene, the effective kinetic barrier for net water intercalation is dominated by the large work required to create and sustain this thermodynamically unfavorable confined state, rather than by the diffusion barrier along the interface. We also find that increasing the number of intercalated water molecules gradually reduces this thermodynamic penalty, likely due to the formation of stabilizing hydrogen-bond networks within the confined region. These findings imply that, during vSFG measurements where the samples are immersed in water, water molecules are likely to be intercalated under monolayer graphene to minimize the system’s free energy (Fig. 6b). However, for thicker graphene (n ≥ 4), water intercalation is energetically disfavored and thus unlikely to occur (Fig. 6c). This provides a plausible explanation for the experimentally observed hydrophilic vSFG features of monolayer graphene on a CaF2 substrate. Rather than resulting from wetting transparency, these spectral signatures are more likely attributable to the presence of intercalated water, which modifies the interfacial structure and vibrational response, yielding features typically associated with hydrophilic surfaces.

Finally, we note that the present simulations treat nuclear motion classically. Nuclear quantum effects (NQEs) are known to induce a red shift in the O–H stretching frequency without substantially modifying the overall spectral profile,[119] and could, in principle, influence the free energy difference between intercalated and slab water through changes in zero-point energy under confinement.[120–123] However, the computed free energy difference between intercalated and slab water exceeds 2.0 eV per water molecule. Established studies indicate that NQEs typically alter water binding free energies by only a fraction of the hydrogen bond strength, on the order of several to tens of millielectronvolts.[120, 123] Given that the observed free energy difference is two orders of magnitude larger, it is unlikely that NQEs would qualitatively change the conclusions drawn here. Nevertheless, a quantitative assessment of NQEs in the thermodynamics of water intercalation remains a computational challenge and represents an important direction for future work.

Methods
4.1 Reference data

To parametrize MLIPs using the ACE method [89, 90], one needs to perform DFT calculations of total energies and atomic forces for a sufficiently large number of representative reference configurations. These reference structures were generated using the on-the-fly machine-learning force field (MLFF) method[91–93] as implemented in VASP(version:6.4.0).[128–131] A kernel-based MLIP potential was constructed using a basis set derived from reference structures, which was dynamically expanded during MD simulations through an active learning approach. When the estimated Bayesian error of the predicted forces at a particular MD step exceeded a predefined threshold, the predicted forces were considered inaccurate for the current configuration. In such cases, the energy and forces of the current structure were recalculated using DFT to ensure accuracy. The DFT calculated energy and forces were used to expand the basis set of reference data and update the MLIP. The error threshold itself was adaptively updated during the MD simulation based on the average Bayesian error from the preceding ten MD steps. This approach allows the MLIP to efficiently sample the configurational phase space and identify representative structures more effectively than conventional AIMD simulations. We trained new potentials using the ACE method rather than relying on the existing on-the-fly MLSFF potential because the ACE approach offers both higher accuracy and significantly improved computational efficiency.[132] In our benchmark tests for graphene-water systems of the same cell size, MD simulations using the ACE potential were approximately two orders of magnitude faster than those using the on-the-fly MLFF potential. This substantial speed advantage is crucial for performing large-scale and long-timescale MD simulations in the present study.

The interaction between ions and electrons was modelled with the projected augmented wave method[133] with an plane-wave energy cutoff of 600 eV. The convergence criterion for electronic self-consistent steps was set to 10−5 eV. Exchange-correlation effects were treated using the Becke-Lee-Yang-Parr (BLYP) functional,[134, 135] combined with D4 dispersion corrections including Becke-Johnson damping[136]. We employed the BLYP-D4 functional for two main reasons. First, a benchmark study by Brandenburg et al.,[137] systematically evaluated the interaction of a single water molecule with different carbon-based materials using various functionals. Among all the methods considered, BLYP-D4 (GGA), TPSS-D4 (meta-GGA),[138] rev-vdW-DF2 (vdW-DF),[139] and PBE0-D4 (hybrid DF)[140] were identified as the most accurate within their respective functional classes. Considering the large size of our simulation cell (at least hundreds of atoms), BLYP-D4 represents the most computationally feasible choice for our simulations, whereas meta-GGA and hybrid functionals are known to be significantly more computationally demanding. Second, Ohto et al.[141] compared different functionals (BLYP-D3, PBE-D3,[142, 143] revPBE-D3,[144] and BLYP) for reproducing vSFG spectra of the water/air interface and found that PBE-D3 and revPBE-D3 yielded a pronounced red shift in the dangling O–H stretch. In contrast, BLYP provided a more reasonable frequency even without accounting for nuclear quantum effects. This suggests that BLYP’s improved performance originates from a partial error cancellation between the missing nuclear quantum effects and the self-interaction correction error.[119] Brillouin zone integrations were performed using a k-point mesh corresponding to a reciprocal-space resolution of 0.25 Å−1. All on-the-fly MLFF MD simulations were performed in the NVT ensemble at a target temperature of 300 K, controlled by the Langevin thermostat.[145] Although a 1.0 fs timestep, close to the upper limit for systems containing O–H stretching vibrations, was used in these short MD simulations, the resulting trajectories were employed only for training data generation and were not utilized for interfacial structure analysis or vSFG spectral simulations.

4.2 Parameterization of potential

The energies, forces, and coordinates of structures calculated with DFT were extracted from the on-the-fly MLFF MD trajectories. We employed the ACE formalism as implemented in the PACEMAKER package (version: 0.2.7)[89, 132, 146] because it provides a systematically improvable, symmetry- adapted representation of short-range many-body interactions that is well suited to multi-component systems such as the present C/F/Ca/O/H interfaces. The linear ACE form allows us to control the balance between basis size, maximum body order, and accuracy, and to integrate the potentials seamlessly with an active-learning workflow and efficient large-scale MD in LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator, version: 2 Aug 2023). Five percent of the initial dataset was set aside as the test set, with the remaining data used for training. For the simplest systems L1 and L5, the ACE potential was constructed using 3000 basis functions and 8376 parameters with a maximum body order ν = 6, corresponding to seven-body interactions. For the system S, the ACE potential included 8000 basis functions and 20128 parameters with a maximum body order ν = 6, corresponding to seven-body interactions. For all systems containing substrate, graphene and water (SnLm), the ACE potential consists of 9895 basis functions and 25430 parameters with a maximum body order ν = 5, corresponding to six-body interactions. A cutoff distance of 7.0 Å was applied across all the systems. During ACE learning, a uniform weighting scheme was applied to all the reference structures in the loss function. The initial parameterization of the potential was further refined through an active learning strategy, allowing for systematic improvement of the potential by incorporating high-error configurations encountered during MD simulations.

4.3 Active learning

To sample a broader range of atomic environments, we performed MD simulations in the NVT ensemble at a target temperature of 300 K, employing the canonical sampling through velocity rescaling (CSVR) thermostat[147] to control the temperature. The integration time step was set to 0.5 fs. The extrapolation grade γ of each atom at every 100 simulation steps was calculated using the D-optimality criterion[148]. Structures with γ >5 were selected as extrapolative structures. The MD simulations were performed with the LAMMPS software[149, 150] and the extrapolation grades were calculated with the ML-PACE package[132]. To identify the most representative atomic environments, the MaxVol algorithm[148, 151] was used to construct a new active set from the structures generated in the previous steps. Structures that had multiple atoms entered the new active set were given high priority. The top candidates were selected for single-point energy calculations with DFT(BLYP-D4). The energy, forces, and coordinates of these structures were added to the previous training dataset to retrain the ACE potential. The process of MD simulation, configuration selection, single point energy calculation, and potential retraining was iteratively repeated until no structures exceeded the extrapolation grade cutoff (γ ≤ 5) within a 5 ns of MD steps (Fig. 7a). The details of the active learning steps and the coverage of intercalated water molecules are summarized in Table 2 and Table 3, respectively. The final ACE potential, trained with the above active learning strategy, was subsequently used for a large-system MD simulation.

4.4 MD simulations

MD simulations were performed using the LAMMPS package in combination with the ML-PACE package. Each system was first equilibrated in the NVT ensemble for 2 ns, followed by a 2 ns production run in the NVE ensemble. During equilibration, the CSVR thermostat was used to make the system reach the target temperature of 300 K. A fixed integration time step of 0.5 fs was used throughout both equilibration and production runs. Atomic coordinates were saved every 1 fs for subsequent analysis. The simulation cell parameters and atom counts for each system are summarized in Suppl. Table 4. For the L5 system, carbon atoms in the bottom two layers were kept fixed, while all other atoms were allowed to move freely. In all the systems involving the CaF2 substrate, the atoms in the bottom layer of the CaF2 slab were frozen to mimic bulk-like conditions, whereas all remaining atoms were fully flexible during the simulations. Rashmi et al.[57] investigated the graphene/water interface using MD simulations with and without including NQEs. Their results showed that NQEs do not significantly alter the structural or hydrogen-bonding topologies of interfacial water molecules, but primarily induce a systematic red shift in the O–H stretching band. Because NQEs mainly leads to a red shift without changing the overall spectral profile, and because the BLYP functional already yields a reasonable O–H stretching frequency owing to error cancellation between the missing NQEs and self-interaction errors,[119] we did not explicitly include NQEs in our MD simulations.

4.5 vSFG spectra simulations

Because the direct contributions of the CaF2 substrate and graphene to vSFG spectra are negligible,[84, 152–154] the experimentally measured vSFG signals originate predominantly from interfacial water molecules. Consequently, the reliability of simulated vSFG spectra mainly depends on the precise evaluation of the dipole moments and polarizability tensors of individual water molecules. The many-body dipole moment surface (MB-µ) and polarizability surface (MB-α) have been demonstrated to reproduce the infrared and Raman spectra of liquid water accurately [155–157], as well as the vSFG spectra of the water/air interface [57, 99, 158]. In this framework, any total property p of an N-molecule system can be decomposed into a many-body expansion, where contributions arise from one-body (1B), two-body (2B), and higher-order interaction terms[159]

pN-mer(1,2,···,N)= ∑i p1B(i)+ ∑i<j p2B(i,j)+··· (1)

where p1B, p2B, ···, refers to the one-body, two-body, ···, contributions, respectively. The 1B term depends solely on the nuclear coordinates of an individual monomer, indexed by i. The 2B contribution involves the nuclear coordinates of the dimer formed by monomer i and j. Therefore, the total dipole moment of the molecular system can be expressed as

µN-mer(1,2,···,N)= ∑i µ1B(i)+ ∑i<j µ2B(i,j)+··· (2)

where µnB denotes the nth interaction-induced dipole term. A similar expansion can be formulated for the system’s polarizability

αN-mer(1,2,···,N)= ∑i α1B(i)+ ∑i<j α2B(i,j)+··· (3)

where αnB represents the nth interaction-induced polarizability. The 1B dipole and polarizability are defined as the corresponding monomer quantities

µ1B =µ(i) (4)
α1B =α(i) (5)

All higher-order terms are defined recursively using Equations 2 or 3. For instance, the 2B dipole term is obtained by subtracting the 1B contributions from the total dipole moment of the dimer

µ2B =µ(1,2)− ∑i=1 2 µ1B(i) (6)

Based on Equations 2 and 3, researchers in Francesco Paesani’s group developed the many-body representations of the dipole moment (MB-µ) and polarizability (MB-α) surfaces for water, derived from highly correlated electronic structure calculations. Detailed descriptions of the training and validation of the MB-µ and MB-α potentials can be found in the work of Medders et al ..[160]

Using our trained ACE potential, which inherently includes the influence of graphene and the substrate’s polarizable effects on the orientation of interfacial water molecules, we first perform MD simulations in the NVE ensemble to generate the required trajectories. The LAMMPS code, in combination with the MBX package (version: 1.2.0), reads the coordinates of all water molecules and computes both their permanent and induced dipole moments using the MB-µ potential. Based on these coordinates and dipole moments, the polarizability of each water molecule is subsequently evaluated using the MB-α potential. For the dipole moments, only the one-body and N-body effects are considered and one-body for polarizability. We acknowledge that the MB-µ and MB-α do not explicitly incorporate the induced dipole moments or atomic polarizabilities in water molecules that arise from interactions with graphene or substrate atoms. Instead, the dipole moment and polarizability of each water molecule depend solely on the geometry and orientation of water molecules obtained from the trajectories generated by our ACE potential. This strategy—considering only the geometry and orientation of water molecules for vSFG spectra simulations—has been widely adopted by other research groups studying solid/liquid interfaces and shows good agreement with the experimental measurements.[55, 57, 68, 86, 161–163]

Once the dipole moments and polarizabilities of all water molecules are obtained, the total time correlation function (TCF) can be written as the sum of the autocorrelation term for the ith water molecule and the cross-correlation term between ith and jth water molecules, respectively.

⟨αbc(t)µa(0)⟩= ∑i ⟨αi,bc(t)µi,a(0)⟩ + ∑i≠j ⟨αi,bc(t)µj,a(0)⟩ (7)

where αbc and µa are the polarizability and dipole moment components of each water molecule, respectively. ⟨···⟩ represents the thermal average. To accurately evaluate the cross-correlation contributions within the limited length of the MD trajectory, a truncated TCF was employed.[164] A cutoff rt was introduced that only water molecule pairs (i and j) separated by less than rt were included in the cross-correlation term. Here, the cutoff radius rt was set to 4.0 Å.

⟨αbc(t)µa(0)⟩= ∑i ⟨αi,bc(t)µi,a(0)⟩ + ∑i≠j gt(rij;rt)⟨αi,bc(t)µj,a(0)⟩ (8)

gt(rij;rt)= {1, rij ≤rt
0, rij >rt (9)

For the slab model, it has two interfaces (water/air and water/graphene or CaF2). Since the two interfaces in the simulation cell are orientated in opposite directions, and only the water/graphene or water/CaF2 interface is of interest, the analysis is restricted to that specific interfacial region. To minimize the influence of the water/air interface on the graphene/water or water/CaF2 interface, a sufficient large number of water molecules was included to create a bulk-like region separating these two interfaces. A screening function gsc(zi)[119, 164] was introduced to the truncated TCF to exclude the contribution from the water/air interface:

⟨αbc(t)µa(0)⟩= ∑i ⟨αi,bc(t)µi,a(0)⟩gsc3(zi)
+ ∑i≠j gt(rij;rt)⟨αi,bc(t)µj,a(0)⟩gsc3(zi) (10)

gsc3(zi)= {1, zi ≥z1
cos2 [π(zi−z2) / 2(z1−z2)], z2 ≤zi <z1
0, zi <z2 (11)

where zi is z coordinate of the ith water molecule. The origin point in the z direction is set to the center mass of the water system(Fig. 7b). In the region zi ≥ z1, the contributions from the dipole moment and polarizability are fully included in the truncated TCF. In the transition region z2 ≤ zi < z1, these contributions are smoothly switched off. For other region, the contributions are set to zero.

The second-order susceptibility χ(2) can be written as a Fourier transform of the truncated TCF of the dipole moment and polarizability components,

χ(2)cba(ω)=i/kBT ∫0 ∞ dteiωt⟨αbc(t)µa(0)⟩ (12)

where ω is the vibrational frequency, i is the imaginary unit, kB is Boltzmann’s constant, a, b, and c are the polarization components of the infrared IR, visible, and sum frequency beams, respectively. Since the vSFG spectra were measured with SSP polarization combination,[35] for the simulation, we only consider the XXZ polarization combination.[45] The imaginary (Im(χ(2))) and real (Re(χ(2))) part of the vSFG spectra were computed by calculating the cosine and sine components of the squared spectrum at different surface depths.

4.6 Free energy calculations

The Gibbs free energy difference (∆G) is given by:

∆G=∆H −T∆S =∆U +p∆V −T∆S ≈∆U −T∆S (13)

where H is the enthalpy, T is the temperature, S is the entropy, U is the internal energy, p is the pressure, and V is the volume. Since the volume change (∆V) is negligible in our system, the Gibbs free energy difference can be approximated as ∆G≈∆U−T∆S. The internal energy is extracted directly from LAMMPS MD simulations. To calculate the entropy of intercalated and slab water molecules, we used the two-phase thermodynamic (2PT) model[165] implemented in DoSPT software (version 0.2.2)[166]. The internal energy and temperature used for free energy calculations correspond to the average values obtained over the first nanosecond of the 2 ns production simulation (Table 8). The trajectory from the first nanosecond was also used as an input for the entropy calculations with DoSPT.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This work combines first-principles molecular dynamics with Atomic Cluster Expansion (ACE) machine learning potentials to systematically investigate the wetting behavior at graphene–water interfaces. The authors construct a range of model systems that include different numbers of graphene layers, a CaF2 substrate, and intercalated water molecules. Through simulations of vibrational sum-frequency generation (vSFG) spectra and free-energy calculations, they explore in depth the effects of the substrate, layer-number dependence, and water intercalation on interfacial properties. The study concludes that graphene is intrinsically hydrophobic and does not exhibit wetting transparency. The experimentally observed hydrophilic behavior of monolayer graphene on hydrophilic substrates is primarily attributed to the combined contribution of water on both sides of the graphene, especially the spontaneous intercalation of water between monolayer graphene and the hydrophilic CaF2 substrate. The results confirm that pristine graphene is hydrophobic and that hydrophobicity becomes more pronounced with increasing graphene thickness. While monolayer graphene readily undergoes water intercalation when prepared under ambient humidity or immersed in water during vSFG measurements, this process becomes energetically unfavorable for multilayer graphene, making intercalation unlikely under the same conditions. However, the study has several issues that merit attention:\n1. The authors adopt a strategy of training separate ACE potentials for different systems. It is recommended to clearly state in the main text or Discussion the motivation and quantitative justification for using multiple system-specific potentials instead of a single “universal” potential. While the current approach may ensure high accuracy within each subsystem, the lack of cross-system transferability and robustness checks could limit the comparability and generalizability of conclusions across systems. It is suggested to add cross-system validation results.\n2. The training data are not described in sufficient detail, making it difficult for readers to assess whether the models adequately sample the relevant atomic environments, and impacting reproducibility and verifiability. It is necessary to provide the relevant datasets and training files in the Supporting Information.\n3. The study investigates how intercalated water between the substrate and graphene varies with graphene layer number and coverage, but does not consider possible interlayer water intercalation between graphene sheets on the water-facing side, nor does it provide a clear rationale for neglecting this scenario. It is need to clarify in the main text or Discussion the physical basis and conditions under which this assumption holds (e.g., in defect-free, edge-sealed multilayer graphene, the interlayer spacing and impermeability make water ingress difficult, requiring edge or defect pathways), and provide supporting analysis or literature evidence. This will clarify the model boundaries and help readers understand the scope of applicability of the conclusions.\n4. The manuscript cites reference 72 (currently listed as Zhao et al., 2014) for key values regarding the migration barrier (~0.19 eV) and diffusion coefficient of water under monolayer graphene. Upon checking, reference 72 appears to discuss ferroelectric ice phases and does not contain these specific kinetic parameters. The authors must carefully verify and correct this citation. The correct reference(s) for these values should be provided, as they are crucial for supporting the argument of facile water intercalation. Furthermore, following the authors' conclusion that water intercalation becomes thermodynamically unfavorable in four-layer graphene, it would be insightful to comment on the expected kinetic barrier in this case.\n5. The authors model the CaF2 substrate using its (111) surface without clarification of terminated-atom. Besides, it is known that different exposed facets (e.g., (100), (110)) can have different surface energies and hydrophilicities. Could the authors please briefly comment on why the (111) surface was chosen and whether the conclusions are expected to hold for other prevalent facets of CaF2?\n6. The strategy of training multiple system-specific potentials is practical for managing complexity. However, it would be valuable to briefly discuss the transferability and expected limitations of these potentials. For instance, if a potential trained on a system with a specific number of intercalated water molecules were applied to simulate a structure with a significantly different water coverage, what would be the expected impact on accuracy?\n7. The authors successfully employed the Atomic Cluster Expansion (ACE) method to develop accurate MLIPs for complex interfaces. It would be helpful for a broader audience if the authors could briefly comment on their specific rationale for choosing ACE over other established MLIP frameworks. What are the advantages of ACE for this particular system, such as its efficiency in describing diverse chemical environments (C, F, Ca, O, H) or its performance for property reconstruction.\n8. The authors highlight that their work advances beyond Ref. 84 by explicitly modeling water on both sides of the graphene sheet—the intercalated layer and the overlying water slab (comprising 512 molecules). This dual-interface approach is a more realistic representation of experimental conditions and is a strength of the study. However, it is not demonstrated that the thickness of this water slab is sufficient to ensure that the structure and dynamics of the water-graphene interface are converged. Besides, the Introduction states that vSFG is surface-specific because the bulk water is centrosymmetric and thus does not generate a signal. This seems to raise a conceptual question: if the bulk is silent, what is the specific physical role of including such a thick water slab in the simulation?\n9. Current SHG simulations require reading static self-consistent wavefunctions to extract band information. For large systems, such calculations are highly computationally demanding. How did the authors achieve the vSFG calculations?", "reviewer_2": "The manuscript addresses whether graphene’s wettability is intrinsically hydrophobic or “wetting transparent” (i.e., reflecting the substrate’s properties). Using molecular dynamics accelerated by a machine-learning atomic cluster expansion (ML-ACE) potential, the authors simulate water at graphene interfaces, considering monolayer and multilayer graphene on a hydrophilic CaF2 substrate. They compute vibrational sum-frequency generation (vSFG) spectra and free-energy profiles for water intercalation. The study finds that pristine graphene is intrinsically hydrophobic, with hydrophobicity increasing as the number of graphene layers grows.\n\nA key result is that monolayer graphene on a hydrophilic substrate can appear hydrophilic, not because graphene is transparent to wetting, but because a thin layer of water intercalates at the graphene–substrate interface and cancels the hydrophobic spectral signature. On this basis, the authors argue that “wetting transparency” is a misleading concept: graphene itself remains hydrophobic, while the substrate influences measurements indirectly via intercalated water.\n\nBelow are some comments that it would be useful to consider before publication:\n1) In a recent study [ACS Nano 19, 4876 (2025)], Rashmi et al. performed advanced simulations of graphene/water using the MB-pol potential with path-integral molecular dynamics, examining interfacial structure and the vSFG spectrum at a neutral graphene surface. Their simulated vSFG spectrum agreed closely with the experimental spectrum of Wang et al. [Angew. Chem. Int. Ed. 63, e202319503 (2024)], supporting the hydrophobic nature of graphene. The authors should discuss their monolayer graphene results in the context of the analyses reported in that study.\n2) The influence of the substrate on the vSFG spectrum at a graphene/water interface has been addressed experimentally by Wang et al. [Nature 615, E1 (2023), Ref. 83]. While the analyses presented in this manuscript are careful and detailed, the conclusions confirm established understanding in the solid/liquid interfaces community. The authors should sharpen the novelty claim and make explicit what goes beyond prior experimental and theoretical work.\n3) The ML potential was trained on BLYP-D4. BLYP is a dated GGA functional with known deficiencies for liquid water (see Gillan et al., J. Chem. Phys. 144, 130901 (2016)). The manuscript should acknowledge these limitations and discuss sensitivity to the choice of functional. If feasible, the authors should bolster confidence by comparing a few key observables (e.g., interfacial orientations, surface tension, contact angle) using a more accurate model. At minimum, the authors should explain why BLYP-D4 was chosen over more modern functionals and summarize its known limitations for water.\n4) All simulations treat nuclei classically. Nuclear quantum effects (NQEs) can affect interfacial structure and the OH-stretch spectrum, and prior path-integral studies have shown broadening/softening and changes in hydrogen-bonding [e.g., Chem. Rev. 116, 7529 (2016)]. Since the authors attribute some deviations from experiment to neglected NQEs, conclusions should be tempered accordingly. A discussion on whether including NQEs would plausibly enhance or reduce intercalation under graphene and how this might impact the main conclusions would also be useful.\n5) Using an ML potential to access larger length/time scales is sensible, but more quantitative validation is needed. Beyond the SI correlation plots for energies/forces, the authors should provide validation on relevant water properties (e.g., density, RDFs, surface tension, and contact angle). This will help establish that the ML-ACE potential faithfully reproduces the reference physics for interfacial water and supports the reported analyses.", "reviewer_3": "The authors employed a combined first-principles and machine learning molecular dynamics simulation, and simulated water/multilayer graphene sheets on CaF2, whose water structure has been debated for a long time. The simulation results showed that “wetting transparency” of the graphene should be explained by the intercalated water between graphene and CaF2, and highlighted the hydrophobic nature of graphene, which is a very important suggestion. Their simulations are thoroughly done and the manuscript is well written. Therefore, I believe that the manuscript can contribute to the understanding of the structure and dynamics of interfacial water, and has a potential interest for the broad readership. However, some parts of the manuscript should be revised to meet the high standards of Nature Communications.\n・Are the MB-mu and MB-alpha, which may be applicable for the water/air interface, applicable for water/solid interfaces? I mean, do the dipole moment surface and polarizability surface cover the water structure modulated by the solid surface? Such a validation would be required.\n・The SFG spectrum for the water/air interface looks a bit deviated from a spectrum computed similarly (J. Phys. Chem. B 122, 4356 (2018).). The presented spectrum (Fig. 2c), in particular, the feature of the 3600 shoulder peak, looks like the low-temperature spectrum (238K) in the reference. Is it because the MLIP is learned from water/graphene interface?\n・If the authors would like to suggest that the experimentally observed vSFG signal arose from signal cancellation, it would be better to compare computed signals with some measured signals and discuss.\n・The timestep of the DFT-MD simulation (1.0 fs) is a bit aggressive for systems containing O-H stretch modes. It is better to show total energies of DFT-MD in addition to those of MLIP-MD. Also, how long DFT-MD simulations are performed is not clearly mentioned.\n・I think that the effect of the cross-correlation term on the spectra is important and should be presented in the supplementary materials.\n・I noticed that the citation (ref. 72) is not consistent.\n・Some labels overlap in Fig. 6a."}, "author_response": "We would like to express our sincere gratude to the three reviewers for their me and thoughul, construcve comments. We have thoroughly revised the manuscript in accordance with their construcve suggesons and have carefully addressed every concern raised. In this document, we provide comprehensive, point-by-point responses to each comment. For clarity, all revisions made to the manuscript have been highlighted with a colored background.\nReviewer #1\n1. The authors adopt a strategy of training separate ACE potenals for different systems.\nIt is recommended to clearly state in the main text or Discussion the movaon and quantave jusficaon for using mulple system-specific potenals instead of a single “universal” potenal. While the current approach may ensure high accuracy within each subsystem, the lack of cross-system transferability and robustness checks could limit the comparability and generalizability of conclusions across systems. It is suggested to add cross-system validaon results.\nWe appreciate the reviewer’s comment. In principle, we fully agree that a single “universal” machine learning interatomic potenal (MLIP) would be highly desirable. In pracce, however, within the current Atomic Cluster Expansion (ACE) framework, this is very challenging for the present systems. ACE is a local, short-ranged expansion: the atomic energy depends only on neighbors within a finite radial cutoff (7 Å in our models) and on a body-order expansion. Consequently, long-range interacons such as electrostacs and dispersion tails are not explicitly included.\nFor the graphene–water systems considered here, the relave energies of different water orientaons at the interface are typically below 0.043 eV (≈ 1 kcal/mol). Accurately capturing these very small energy differences, which control the vSFG response, requires “sub-chemical accuracy” of the potenal. Mixing all environments—bare CaF₂, suspended graphene, supported graphene with and without intercalated water, and different graphene layer numbers—into a single ACE model would require a dramacally larger basis and/or body order to reach the same level of precision. In parcular, the relevant water–water separaons differ strongly between S L and S L systems: intercalated and n 1 n 4 slab waters fall inside the cutoff in S L , but are separated by 15 Å in S L , i.e., far beyond n 1 n 4 the range of the short-ranged expansion. A universal short-range potenal would ≳ therefore have to interpolate between physically disnct regimes, which we expect to degrade its accuracy for both.\nFor this reason, we opted for system-specific ACE potentials (ACE–L , ACE–L , ACE–S, ACE– 1 5 SL , ACE–SL , ACE–SL , ACE–S L and ACE–S L ), as now stated more explicitly in Sec. 2.2 1 4 5 n 1 n 4 of the revised manuscript (Page 7). Each potential is trained on configurations representative of the corresponding system, while all share common building blocks (bulk liquid water) and are fitted to the same DFT(BLYP-D4) reference. As summarized in Table S3, all models achieve very small energy and force errors (≤ 0.50 meV/atom and < 40 meV/Å on both training and test sets), and long NVE simulations show stable temperatures and total energies (Figs. S11–S29), confirming that the learned potentials are robust for their target systems.\nWe emphasize that the comparisons across systems in the manuscript (e.g., trends in interfacial water orientation and vSFG spectra between L /L and SL /SL /SL , and 1 5 1 4 5 between systems with different intercalated water coverages) are therefore comparisons between DFT-consistent models, rather than extrapolations of a single potential far beyond its training domain. The system-specific potentials are not intended to be transferable across distinct systems. For example, applying ACE-S L to the S L system n 1 n 4 causes the MD simulation to terminate within a few steps, as the extrapolation grade becomes prohibitively large, indicating that the configuration lies far outside the training domain and that the potential cannot guarantee reliable predictions. Extending the present approach toward a more transferable, long-range-corrected ACE potential that covers all these regimes is an important direction for future work but lies beyond the scope of the current study.\nA detailed discussion on the raonale for using system-specific rather than universal potenal has been added to the revised manuscript in secon 2.2. The revisions have been made on pages 7-8 of the manuscript and are highlighted with a colored background.\n2. The training data are not described in sufficient detail, making it difficult for readers to assess whether the models adequately sample the relevant atomic environments, and impacng reproducibility and verifiability. It is necessary to provide the relevant datasets and training files in the Supporng Informaon.\nWe thank the reviewer for emphasizing the importance of a transparent description of the training data and of providing the underlying datasets. Based on the reviewer’s comments, all training and test datasets, the final machine learning interatomic potentials as well as structural files and input files for MD simulations are available on Zenodo (htps://doi.org/10.5281/zenodo.17836612).\nIn the revised SI, we now provide a more detailed description of the training and test set compositions. Specifically, Section 1 of Table S1 summarizes, for each ACE model, the total number of reference structures, the number of configurations originating from the on-the-fly machine-learning force field molecular dynamics (OFMLFFMD) stage and from the active learning molecular dynamics (ALMD) stage, and the lengths of the underlying MD simulations. For the intercalated water systems, the range of water coverages sampled during OFMLFFMD and ALMD is listed in Table S2, covering the conditions studied in the production simulations. Table S3 reports, for each system, the size of the training and test sets together with the energy and force RMSEs; each training set includes 1843 bulk-water supercells to ensure a consistent description of liquid water across all systems. The revisions have been made on pages 3-4 of the SI.\n3. The study invesgates how intercalated water between the substrate and graphene varies with graphene layer number and coverage, but does not consider possible interlayer water intercalaon between graphene sheets on the water-facing side, nor does it provide a clear raonale for neglecng this scenario. It is needed to clarify in the main text or Discussion the physical basis and condions under which this assumpon holds (e.g., in defect-free, edge-sealed mullayer graphene, the interlayer spacing and impermeability make water ingress difficult, requiring edge or defect pathways), and provide supporng analysis or literature evidence. This will clarify the model boundaries and help readers understand the scope of applicability of the conclusions.\nThank you for this comment. This theorecal work primarily aims to interpret the experimentally measured vSFG spectra and to elucidate the disnct variaons in vibraonal signatures observed with increasing numbers of graphene layers, as reported by Kim et al..1 In their study, the defect density (defects per carbon atom) for graphene with different layer numbers was determined to be 13.6±3.8 ꓫ 10-6 for monolayer graphene, 9.3±2.6 ꓫ 10-6 for four-layer graphene, and 5.1±1.5 ꓫ 10-6 for five-layer graphene. Given these extremely low defect densies, we model the graphene sheets as prisne and defect-free, focusing on the interior region of large-area flakes where edge effects are negligible.\nIn this regime, it is well established that prisne graphene is effecvely impermeable to water molecules.2 Consequently, for mullayer prisne graphene, water cannot intercalate between the graphene layers from the water-facing side. Arif et al.3 used atomic force microscopy (AFM) to measure the thickness of mullayer graphene under relave humidies ranging from 10% to 80% and observed only minimal changes in thickness, concluding that humidity-induced water intercalaon into mullayer graphene is negligible. Based on these experimental findings, water intercalaon between graphene layers was not considered in our simulaons.\nOur simulaon cells employ periodic boundary condions and contain no explicit edges; they thus represent the interior, defect-free, edge-sealed poron of mullayer graphene. Within this framework, water is allowed to intercalate only between the botom graphene layer and the hydrophilic CaF substrate, where capillary forces and the favorable free energy we compute make intercalaon thermodynamically likely for the monolayer graphene system. In contrast, intercalaon between graphene layers on the water-facing side is excluded, because under the defect-free, edge-sealed condions we consider, such a process is both kinecally and structurally unlikely on the experimental mescale.\nWe have revised the manuscript to make this assumpon and its physical basis explicit, and to clarify that our conclusions regarding the absence of weng transparency in mullayer graphene apply to high-quality, large-area films with very low defect densies and effecvely sealed edges. Possible effects of interlayer water in more defecve or edge-dominated samples are outside the scope of the present work and are now explicitly noted as such in the Discussion. The revisions have been made on page 13 and page 19 of the manuscript and are highlighted with a colored background.\n4. The manuscript cites reference 72 (currently listed as Zhao et al., 2014) for key values regarding the migraon barrier (~0.19 eV) and diffusion coefficient of water under monolayer graphene. Upon checking, reference 72 appears to discuss ferroelectric ice phases and does not contain these specific kinec parameters. The authors must carefully verify and correct this citaon. The correct reference(s) for these values should be provided, as they are crucial for supporng the argument of facile water intercalaon. Furthermore, following the authors' conclusion that water intercalaon becomes thermodynamically unfavorable in four-layer graphene, it would be insighul to comment on the expected kinec barrier in this case.\nWe sincerely thank the reviewer for the careful observaon. We have corrected the inconsistency in the citaon of this reference 72, which is now listed as reference 70 (Front. Phys. 9:740627). All other references have been carefully reexamined to ensure accuracy and consistency throughout the manuscript.\nIn the referenced paper, the author calculated the migraon barrier for a single water molecule moving from the top of a Ca atom to a hollow site beneath monolayer graphene. The forward and backward barriers were reported to be 0.19 eV and 0.25 eV, respecvely. They also obtained a high diffusion coefficient for water molecules at room temperature (300 K), 8.35 ꓫ 10-6 cm2/s, indicang rapid diffusion. This is one of the reasons that movated us to consider the possibility of water intercalaon between monolayer graphene and the substrate. Water intercalaon between monolayer graphene and hydrophilic substrates has been extensively studied,4-10 although only a few studies have invesgated the corresponding vSFG spectra of intercalated water molecules.11 Previous works have clearly demonstrated that water can spontaneously intercalate between monolayer graphene and hydrophilic substrates.8, 11\nIn our free-energy calculaons for the intercalated water systems, the average free energy of water in the water slab of each system was set to zero. We then compared the average free energy of the interfacial water molecules to that reference. The results show that water intercalaon is thermodynamically more favorable than adsorpon on top of monolayer graphene, whereas the opposite trend is observed for the four-layer graphene system. This behavior remains consistent across different interfacial coverages. We did not perform kinec barrier calculaons for the intercalaon process because, in our model, water molecules are present on both sides of the graphene sheet. In the simplified model, where only a single water molecule beneath the graphene is considered, and no water is present on the opposite side (as in Front. Phys. 9:740627), the interacon between the graphene and slab water is neglected. This is fundamentally different from our present model. Restricng the analysis to water beneath the graphene would not be consistent with the thermodynamic model used in this study.\nWe have replaced the sentence:\n“Furthermore, due to the relavely low migraon barrier (approximately 0.19 eV) and high diffusion coefficient (8.35 × 10⁻⁶ cm²/s) of water,[72] it is highly likely to have intercalated water molecules between the hydrophilic substrate and graphene.”\nwith\n“Furthermore, first-principles calculaons for the CaF₂(111)/graphene system have shown that once a water molecule is confined between monolayer graphene and the CaF₂(111) substrate, its lateral diffusion between neighboring adsorpon sites proceeds with a very low migraon barrier of approximately 0.19 eV, leading to a diffusion coefficient of about 8.35 × 10⁻⁶ cm²/s at 300 K.[70] Such facile diffusion implies that, once intercalated, water molecules can rapidly spread under the graphene layer, forming an extended confined layer.”(Pages 13-14 of the manuscript)\nAer the following sentence\n“In contrast, the free energy of an intercalated water molecule becomes posive (> 1.5 eV) for the four-layer graphene system (Fig. 6a), indicang an energecally unfavorable intercalaon process.”\nwe have added a clarifying sentence:\n“At room temperature, a posive free-energy penalty of about 1.7–2.1 eV per intercalated water molecule (Table S8) corresponds to an exponenally small equilibrium populaon of confined water under four-layer graphene, even though the local lateral diffusion barrier within the CaF₂(111)/graphene gap is expected to remain low. Thus, for mullayer graphene, the effecve kinec barrier for net water intercalaon is dominated by the large work required to create and sustain this thermodynamically unfavorable confined state, rather than by the diffusion barrier along the interface.” (Page 19 of the manuscript)\n5. The authors model the CaF2 substrate using its (111) surface without clarificaon of the terminated-atom. Besides, it is known that different exposed facets (e.g., (100), (110)) can have different surface energies and hydrophilicies. Could the authors please briefly comment on why the (111) surface was chosen and whether the conclusions are expected to hold for other prevalent facets of CaF2?\nThank you for this comment. It is indeed true that different facets of CaF₂ exhibit disnct surface energies and adsorpon energies for water molecules. Cen et al.12 calculated the surface energies of CaF₂(111), CaF₂(110), and CaF₂(210) surfaces using the DFT (PBE funconal), obtaining values of 0.52, 0.82, and 1.11 J/m², respecvely. The adsorpon energies of a single water molecule on CaF₂(111) and CaF₂(110) are -0.50 eV and -0.20 eV, respecvely. For CaF₂(210), there are two possible adsorpon sites for a single water molecule, with adsorpon energies of -0.80 eV and -1.31 eV, respecvely.\nThe CaF₂(111) surface is chemically inert because the surface is fluorine terminated, resulng in stability in air, low solubility in water, and the absence of dangling bonds or - charge traps.13 This inert surface configuraon has been shown to effecvely reduce interfacial interacons between 3D bulk materials and 2D materials, giving rise to a quasi-van der Waals gap.13 The resulng quasi-van der Waals interacon between 2D materials and CaF₂ produces an ultraclean interface, where carrier moon is largely free from impurity scatering and charge traps.14 This explains why the CaF₂(111) surface is commonly used in the fabricaon of van der Waals devices. For these reasons, many studies have adopted the CaF₂(111) surface as the default for both experimental vSFG spectroscopy measurements and simulaons of vSFG spectra for CaF₂/water and CaF₂/graphene/water interfaces.1, 11, 15-17\nSince the surface and adsorpon energies differ between facets, the extent to which different CaF₂ facets influence the vSFG spectra remains uncertain. Nevertheless, provided that the substrate remains strongly hydrophilic, the qualitave conclusions reached here—the intrinsic hydrophobicity and non-weng transparency of prisne graphene, and the dominant role of intercalated water in the apparent hydrophilicity of monolayer graphene—are expected to remain valid. We have added sentences in the manuscript to discuss the potenal influence of different facets. We have clarified the fluorine-terminated (111) surface on page 10 and added further discussion on the effect of facets on page 21. All changes are highlighted with a colored background.\n6. The strategy of training mulple system-specific potenals is praccal for managing complexity. However, it would be valuable to briefly discuss the transferability and expected limitaons of these potenals. For instance, if a potenal trained on a system with a specific number of intercalated water molecules were applied to simulate a structure with a significantly different water coverage, what would be the expected impact on accuracy?\nWe appreciate the reviewer’s request for clarificaon. For the intercalated water systems S L and S L , we trained two separate potenals, ACE-S L and ACE-S L . It is true that n 1 n 4 n 1 n 4 ACE-S L can only be used for molecular dynamics simulaons of the S L system and n 1 n 1 cannot be applied to S L , and vice versa for ACE-S L . If ACE-S L is forcibly applied to the n 4 n 4 n 1 S L system, the MD simulaon terminates within a few steps because the extrapolaon n 4 grade becomes extremely large, indicang that the configuraon lies far outside the potenal’s training domain, which cannot guarantee the accuracy of the potenal. The explanaon has been added on pages 7–8 of the manuscript and is highlighted with a colored background.\nTable 1 Number of water molecules per nm2 in SnL1 and SnL4 systems during on-the-fly machine-learning force field molecular dynamics (OFMLFFMD) and acve learning molecular dynamics (ALMD) simulaons.\nS L S L\nn 1 n 4\n1.00, 2.01, 8.04,\nOFMLFFMD 1.00, 4.02, 8.04, 16.08\n12.06, 16.08\n0.50, 1.00, 2.01, 2.51, 3.01,\n0.50, 1.00, 1.51, 4.02, 6.03,\nALMD 4.02, 6.03, 8.04, 10.05,\n8.04, 16.08, 19.09\n12.06, 16.08, 23.11, 36.18\nDuring the training process, even though relavely small cell sizes were used for training data preparaon, we considered a wide range of water coverages to ensure that ACE-S L n 1 and ACE-S L perform well across different coverages within their respecve systems. The n 4 specific water coverage values used in the on-the-fly machine-learning force field molecular dynamics (OFMLFFMD) and acve learning molecular dynamics (ALMD) simulaons are listed in Table 1. This table has been added to SI (Table S2, Page 4 of SI).\nThese coverage values encompass the intercalated water coverage condions adopted in our long-me and large-scale MD simulaons. Therefore, we are confident that the present simulaons at the intercalated water coverage are reliable. All training datasets associated with the various water coverages are available on Zenodo (htps://doi.org/10.5281/zenodo.17836612).\nTo directly answer the reviewer’s “what happens if…” question and to delineate the domain in which the potentials can be trusted, we have added the following sentences in SI:\n“For the intercalated-water systems S L and S L , separate ACE potentials (ACE– n 1 n 4 S L and ACE–S L ) were trained. In both cases, the on-the-fly machine-learning n 1 n 4 force field molecular dynamics (OFMLFFMD) and active learning molecular dynamics (ALMD) stages sampled a broad range of intercalated water coverages (from ~0.5 up to ~36 water molecules per nm² for S L and up to ~19 water n 1 molecules per nm² for S L ), as summarized in Table S2. All production simulations n 4 reported here use coverage within these training ranges; applying the present potentials to substantially higher coverages or qualitatively different confined water structures would constitute extrapolation and may require additional training data.” (Page 4 of the SI)\n7. The authors successfully employed the Atomic Cluster Expansion (ACE) method to develop accurate MLIPs for complex interfaces. It would be helpful for a broader audience if the authors could briefly comment on their specific raonale for choosing ACE over other established MLIP frameworks. What are the advantages of ACE for this parcular system, such as its efficiency in describing diverse chemical environments (C, F, Ca, O, H) or its performance for property reconstrucon.\nThank you for this comment. To provide some background, this project was iniated at the beginning of 2023. At that me, we carefully reviewed the pioneering ACE studies by Ralf Drautz, Yury Lysogorskiy, and Anton Bochkarev.18-20 Their work demonstrated that the ACE method provides a systemacally improvable and symmetry-compliant framework for represenng short-range many-body interacons, offering an excellent balance between accuracy and computaonal efficiency for mul-component systems. The comparave accuracy of the ACE potenal relave to other exisng MLIPs is illustrated in Figure 1. Furthermore, the pretrained ACE potenal can undergo iterave refinement through acve learning, enabling systemac improvement across diverse atomic configuraons. In addion, the ACE potenal is computaonally efficient: a single force evaluaon requires only 0.32 ms/atom and 0.80 ms/atom for the Cu and Si ACE models, respecvely (mings were obtained on a single core of an Intel(R) Xeon(R) Gold 6132 CPU using GCC 7.3.0 and the LAMMPS version from 4 February 2020).18 These speeds are sufficiently fast for large-scale MD simulaons. Overall, the combinaon of high accuracy, computaonal efficiency, and robust acve learning capability movated our decision to employ the ACE method for our system.\nTo clarify our choice of the ACE method, we have added the following sentence:\n“We trained new potenals using the ACE method rather than relying on the exisng on-the-fly MLFF potenal because the ACE approach offers both higher accuracy and significantly improved computaonal efficiency.[125] In our benchmark tests for graphene/water systems of the same cell size, MD simulaons using the ACE potenal were approximately two orders of magnitude faster than those using the on-the-fly MLFF potenal. This substanal speed advantage is crucial for performing large-scale and long-mescale MD simulaons in the present study.”(Page 23 of the manuscript)\nWe have replaced the sentence:\n“They were used to construct a more accurate and faster MLIP with the ACE method, as implemented in the PACEMAKER package.”\nwith\n“We employed the ACE formalism as implemented in the PACEMAKER package[89, 125, 140] because it provides a systemacally improvable, symmetry-adapted representaon of short-range many-body interacons that is well suited to mul-component systems such as the present C/F/Ca/O/H interfaces. The linear ACE form allows us to control the balance between basis size, maximum body order, and accuracy, and to integrate the potenals seamlessly with an acve-learning workflow and efficient large-scale MD in LAMMPS.” (Page 24 of the manuscript)\nFigure 1 Test RMSE vs. computaonal cost for Cu (a) and Si (b) for ACE potenals compared to other potenals. Figure form Ref.18(CC BY 4.0).\n8. The authors highlight that their work advances beyond Ref. 84 by explicitly modeling water on both sides of the graphene sheet—the intercalated layer and the overlying water slab (comprising 512 molecules). This dual-interface approach is a more realisc representaon of experimental condions and is a strength of the study. However, it is not demonstrated that the thickness of this water slab is sufficient to ensure that the structure and dynamics of the water-graphene interface are converged. Besides, the Introducon states that vSFG is surface-specific because the bulk water is centrosymmetric and thus does not generate a signal. This seems to raise a conceptual queson: if the bulk is silent, what is the specific physical role of including such a thick water slab in the simulaon?\nThank you for this comment. Researchers in Ref. 84 (Nat. Commun. 16, 7288 (2025)) employed heterodyne-detected vSFG spectroscopy to study varying amounts of intercalated water molecules between monolayer graphene and the CaF₂ substrate, providing insights into the structure and properes of water under nanoconfinement. 12 They observed that the vSFG spectra depended sensively on the amount (thickness) of intercalated water molecules.\nTo invesgate the weng transparency of graphene, most vSFG spectroscopy measurements for the graphene/CaF₂ system are performed on samples immersed in water, at least covering the graphene layer and part of the substrate. Our simulaon results indicate that water intercalaon between monolayer graphene and the substrate is energecally favorable, resulng in water molecules present on both sides of the graphene. Consistently, our simulaons show that the vSFG spectra varied with the amount of intercalated water. We therefore believe that these findings movate Operando vSFG spectroscopy measurements. To the best of our knowledge, no experimental studies have yet invesgated vSFG spectra for varying amounts of intercalated water simultaneously with water molecules on the opposite side of graphene.\nTo clarify that, we have replaced the sentence:\n“However, their study did not consider water molecules on the opposite side of the graphene. Our simulaons demonstrate that the coexistence of intercalated and adsorbed water layers induces complex spectral interplay, as evident in the Im(χ(2)) spectra of the S L system (Fig. 4c). This dual-interfacial configuraon H 1 scenario is worth invesgang experimentally, as it may beter reflect realisc condions and improve the understanding of graphene–water interfacial behavior.”\nwith\n“Their results provided important informaon about the structure and properes of nanoconfined water, showing that the vSFG spectra are highly sensive to the thickness of the intercalated water layer. In most studies examining the weng transparency of graphene, vSFG measurements on graphene/CaF are performed 2 with samples immersed in water, ensuring coverage of both the graphene layer and part of the substrate. Our simulaons indicate that water intercalaon between monolayer graphene and the substrate is energecally favorable, leading to the presence of water molecules on both sides of the graphene. The coexistence of intercalated and adsorbed water layers in our simulaons gives rise to complex spectral features, as reflected in the Im(χ(2)) spectra of the S L system (Fig. 4c). n 1 These observaons strongly movate future operando vSFG experiments designed to probe systems with varying amounts of intercalated waters together with water layers on the opposite side of graphene—a configuraon that, to the best of our knowledge, has not yet been experimentally explored.” (Pages 20-21 of the manuscript)\nRegarding the number of water molecules in the slab, our simulaon modules consist of both a graphene/water interface and a water/air interface (Figure 2a). To minimize the influence of the water/air interface on the graphene/water interface, a sufficiently large number of water molecules was included to form a bulk-like region separang these two interfaces. Figure 2b presents the density profiles of water along the Z-axis of the SL 1 system. Although the bulk region does not directly contribute to the vSFG response, it is crucial for suppressing mutual interference between the interfaces, thereby requiring a relavely thick water slab. The corresponding water density relave to the graphene surface is provided in the SI (secon Water density profile along the Z axis, Figure S68-S86).\nFigure 2 a, Atomisc model of the SL 1 system. b, Water density profile along the surface normal (Z-axis), obtained from the last 2 ns of the NVE ensemble molecular dynamics trajectory for the SL system. The mean Z-coordinate of the 1 monolayer graphene carbon atoms was defined as the reference plane (z = 0).\n9. Current SHG simulaons require reading stac self-consistent wavefuncons to extract band informaon. For large systems, such calculaons are highly computaonally demanding. How did the authors achieve the vSFG calculaons?\nThe 1B dipole and polarizability are defined as the corresponding monomer quanes 1𝐵𝐵 𝜇𝜇 (𝑖𝑖) = 𝜇𝜇(𝑖𝑖) (4) 1𝐵𝐵 All higher-order terms are defined rec𝛼𝛼ursi ( v𝑖𝑖e ) ly= us𝛼𝛼in (𝑖𝑖g ) Equaons 2 or 3. For instance, the ( 25B ) dipole term is obtained by subtracng the 1B contribuons from the total dipole moment of the dimer 15 2 2𝐵𝐵 1𝐵𝐵 𝜇𝜇 (1,2) = 𝜇𝜇(1,2)−𝜇𝜇 (𝑖𝑖) (6) 𝑖𝑖\nBased on Equaons 2 and 3, researchers in Francesco Paesani’s group developed the many-body representaons of the dipole moment (MB-μ) and polarizability (MB-α) surfaces for water, derived from highly correlated electronic structure calculaons. Detailed descripons of the training and validaon of the MB-μ and MB-α potenals can be found in the work of Medders et al..22\nUsing our trained ACE potenal, we first perform molecular dynamics simulaons in the NVE ensemble to generate the required trajectories. The LAMMPS code, in combinaon with the MBX package, reads the coordinates of all water molecules and computes both their permanent and induced dipole moments using the MB-μ potenal. Based on these coordinates and dipole moments, the polarizability of each water molecule is subsequently evaluated using the MB-α potenal. Once the dipole moments and polarizabilies of all water molecules are obtained, the truncated me correlaon funcon (TCF; Equaon 5 in the manuscript) can be computed. The vSFG spectra are obtained from the Fourier transform of the truncated TCF. This approach has been successfully employed in several previous studies of vSFG spectra.23-25 To further clarify the vSFG calculaon process, we have added the relevant equaons and detailed explanaons to the SI on pages 88–89.\n1) In a recent study [ACS Nano 19, 4876 (2025)], Rashmi et al. performed advanced simulaons of graphene/water using the MB-pol potenal with path-integral molecular dynamics, examining interfacial structure and the vSFG spectrum at a neutral graphene surface. Their simulated vSFG spectrum agreed closely with the experimental spectrum of Wang et al. [Angew. Chem. Int. Ed. 63, e202319503 (2024)], supporng the hydrophobic nature of graphene. The authors should discuss their monolayer graphene results in the context of the analyses reported in that study.\nWe thank the reviewer for drawing our atenon to the recent work of Rashmi et al. (ACS Nano 19, 4876–4886 (2025)). Using the MB-pol many-body potenal combined with path-integral molecular dynamics and MB-µ/MB-α based vSFG simulaons, they showed that the neutral suspended graphene/water interface exhibits a pronounced dangling OH peak in Im χ(2) that closely matches the HD-vSFG spectrum of Wang et al., thereby confirming the intrinsically hydrophobic character of prisne graphene.\nOur monolayer system (L ) corresponds to the same physical situaon of a defect-free, 1 charge-neutral graphene sheet floang on water. As shown in Fig. 2c, the simulated Im χ(2) for L displays (i) a strong posive dangling-OH peak at ~3665 cm⁻¹, and (ii) reduced 1 intensity in the hydrogen-bonded region around 3200–3400 cm⁻¹ relave to bulk-like water, in close analogy to both the experimental spectra of Wang et al. and the quantum MB-pol simulaons of Rashmi et al. This agreement indicates that our ACE-based MD + MB-µ/MB-α framework reproduces the same hydrophobic interfacial structure obtained with the more expensive MB-pol/path-integral approach.\nAt the same me, our work is complementary to Rashmi et al.: while they focused on the dependence of the graphene/water vSFG spectrum on surface charge at a suspended interface, we invesgate how a hydrophilic CaF₂ substrate, graphene layer number, and thermodynamically controlled water intercalaon jointly modulate the apparent wetability and vSFG response. Together, these studies strengthen the picture that prisne suspended graphene is intrinsically hydrophobic, and that the experimentally observed 17 “hydrophilic” signatures for supported monolayers originate from substrate and intercalated-water effects rather than from a change in the intrinsic graphene/water interacon. We have revised the Introducon and the secon “Effect of various graphene layers” (Pages 3 and 10) to explicitly discuss our monolayer results in the context of Rashmi et al.’s analyses.\n2) The influence of the substrate on the vSFG spectrum at a graphene/water interface has been addressed experimentally by Wang et al. [Nature 615, E1 (2023), Ref. 83]. While the analyses presented in this manuscript are careful and detailed, the conclusions confirm established understanding in the solid/liquid interfaces community. The authors should sharpen the novelty claim and make explicit what goes beyond prior experimental and theorecal work.\nWe appreciate the reviewer’s comment. We agree that in the referenced work(Nature 615, E1 (2023))26, the authors invesgated the effects of the substrate and intercalated water molecules under different external electrical potenals. At negave potenals (≤ –1.6 V, Fig. 1d in Ref 26), they observed a negave peak around 3630 cm⁻¹, which was atributed to the OH stretching mode of CaOH species. However, this negave peak disappeared at posive or weakly negave potenals (≥ –1.0 V). In our simulaons, which were performed without applying any external potenal, we found that as the number of intercalated water molecules increases, the posive peak corresponding to the dangling OH gradually transforms into a negave peak (Figure 3a). As shown in Figure 3b, this change originates from the intercalated water molecules themselves, as the amplitude of the negave peak increases with their number. This finding movates our recommendaon to perform Operando vSFG spectroscopy measurements to verify this behavior. The corresponding discussion has been added to on page 16 of the manuscript and is highlighted with a background color.\nTo clarify the novelty of our work beyond previous experimental and theorecal studies, we have added the following two paragraphs to the manuscript.\n“Compared with earlier experimental and theorecal studies—typically limited to either suspended graphene/water interfaces or supported monolayer graphene under fixed condions—our work (i) disentangles the individual roles of the substrate, graphene layer number, and intercalated water in governing the vSFG response; (ii) idenfies cancellaon between intercalated and adsorbed water signals as the microscopic origin of the “hydrophilic” vSFG features observed for supported monolayer graphene on CaF ; and (iii) quanfies how the 2 thermodynamic driving force for water intercalaon varies with graphene thickness, thereby explaining the experimentally reported thickness-dependent transion in apparent wetability.” (Page 5 of the manuscript)\n“Our results not only confirm the intrinsic hydrophobicity of prisne mono- and mullayer graphene and the strong hydrophilicity of CaF (111), but also offer a 2 microscopic reinterpretaon of the vSFG spectra of supported graphene. Specifically, we show that (i) a defect-free CaF (111) substrate by itself does not 2 make graphene “weng transparent”; (ii) the seemingly hydrophilic vSFG features of supported monolayer graphene originate from cancellaon between the 19 spectral contribuons of intercalated and adsorbed water, enabled by thermodynamically favorable water intercalaon; and (iii) this intercalaon becomes increasingly unfavorable as the number of graphene layers increases, thereby explaining the experimentally observed thickness-dependent evoluon of the 3600–3700 cm-1 feature.” (Page 20 of the manuscript)\n3) The ML potenal was trained on BLYP-D4. BLYP is a dated GGA funconal with known deficiencies for liquid water (see Gillan et al., J. Chem. Phys. 144, 130901 (2016)). The manuscript should acknowledge these limitaons and discuss sensivity to the choice of funconal. If feasible, the authors should bolster confidence by comparing a few key observables (e.g., interfacial orientaons, surface tension, contact angle) using a more accurate model. At minimum, the authors should explain why BLYP-D4 was chosen over more modern funconals and summarize its known limitaons for water.\nWe appreciate the reviewer’s comment. We agree that the BLYP funconal has well-known limitaons, such as an overstructured hydrogen-bond network and an overesmaon of the bulk water density. We chose to employ the BLYP-D4 funconal for two main reasons. First, in the benchmark study by Brandenburg et al,27 various funconals were evaluated for a single water molecule interacng with different carbon-based structures. Among all the methods considered, BLYP-D4 (GGA), TPSS-D4 (meta-GGA), rev-vdW-DF2 (vdW-DF), and PBE0-D4 (hybrid DF) were idenfied as the most accurate within their respecve funconal classes. Given the large size of our model system (at least hundreds of atoms), BLYP-D4 represents the most computaonally feasible choice for our simulaons, whereas meta-GGA and hybrid funconals are known to be significantly more computaonally demanding. Second, in the theorecal work by Ohto et al28 , different funconals (BLYP-D3, PBE-D3, revPBE-D3, and BLYP) were tested for reproducing the vSFG spectra of the water/air interface. They found that PBE-D3 and revPBE-D3 produced a significant red shi in the dangling O–H stretch vibraonal modes, whereas BLYP yielded a more reasonable frequency without including nuclear quantum effects. This suggests that the accuracy of BLYP for vibraonal frequency arises from an 20 error cancellaon between the missing nuclear quantum effects and the self-interacon correcon error.29\nBased on these consideraons, we employed the BLYP-D4 funconal to perform the molecular dynamics simulaons used to generate our training data. The good agreement between our simulated vSFG spectra of the floang monolayer graphene/water and CaF₂/water interfaces and previously reported experimental and theorecal vSFG spectra further confirms that the BLYP-D4 funconal is well-suited for our system. The following sentences explaining the reason for using BLYP-D4 funconal have been added to the Methods secon 5.1 and are highlighted with a colored background.\n“We employed the BLYP-D4 funconal for two main reasons. First, a benchmark study by Brandenburg et al.,27 systemacally evaluated the interacon of a single water molecule with different carbon-based materials using various funconals. Among all the methods considered, BLYP-D4 (GGA), TPSS-D4 (meta-GGA), rev-vdW DF2 (vdW-DF), and PBE0-D4 (hybrid DF) were idenfied as the most accurate within their respecve funconal classes. Considering the large size of our simulaon cell, BLYP-D4 represents the most computaonally feasible choice for our simulaons, whereas meta-GGA and hybrid funconals are known to be significantly more computaonally demanding. Second, Ohto et al.28 compared different funconals (BLYP-D3, PBE-D3, revPBE-D3, and BLYP) for reproducing vSFG spectra of the water/air interface and found that PBE-D3 and revPBE-D3 yielded a pronounced red shi in the dangling O–H stretch. In contrast, BLYP provided a more reasonable frequency even without accounng for nuclear quantum effects. This suggests that BLYP’s improved performance originates from a paral error cancellaon between the missing nuclear quantum effects and the self-interacon correcon error.29”(Page 23 of the manuscript)\n4) All simulaons treat nuclei classically. Nuclear quantum effects (NQEs) can affect interfacial structure and the OH-stretch spectrum, and prior path-integral studies have shown broadening/soening and changes in hydrogen-bonding [e.g., Chem. Rev. 116, 21 7529 (2016)]. Since the authors atribute some deviaons from experiment to neglected NQEs, conclusions should be tempered accordingly. A discussion on whether including NQEs would plausibly enhance or reduce intercalaon under graphene and how this might impact the main conclusions would also be useful.\nWe appreciate the reviewer’s comment. We agree that nuclear quantum effects (NQEs) can modify both inter- and intramolecular interacons, thereby influencing the vibraonal spectra and causing a red shi in frequency compared to spectra obtained with classical nuclei. Rashmi et al.25 invesgated the graphene/water interface using molecular dynamics (MD) simulaons with and without including NQEs. Their results showed that NQEs do not significantly alter the structural or hydrogen-bonding topologies of interfacial water molecules, they do lead to a systemac red shi of the OH-stretch band. Since NQEs mainly induce a red shi without changing the overall spectral profile, and because the BLYP funconal already provides a reasonable frequency for the OH stretch mode due to an error cancellaon between missing NQEs and self-interacon correcon errors,29 we did not explicitly include NQEs in our MD simulaons. We have included an explanaon on page 26 of the manuscript, highlighted with a colored background, regarding why NQEs were omited from our MD simulaons.\n“Rashmi et al.[57] invesgated the graphene/water interface using MD simulaons with and without including nuclear quantum effects (NQEs). Their results showed that NQEs do not significantly alter the structural or hydrogen-bonding topologies of interfacial water molecules, but primarily induce a systemac red shi in the O–H stretching band. Because NQEs mainly leads to a red shi without changing the overall spectral profile, and because the BLYP funconal already yields a reasonable O–H stretching frequency owing to error cancellaon between the missing NQEs and self-interacon errors,[138] we did not explicitly include NQEs in our MD simulaons.”\n5) Using an ML potenal to access larger length/me scales is sensible, but more quantave validaon is needed. Beyond the SI correlaon plots for energies/forces, the 22 authors should provide validaon on relevant water properes (e.g., density, RDFs, surface tension, and contact angle). This will help establish that the ML-ACE potenal faithfully reproduces the reference physics for interfacial water and supports the reported analyses.\nThank you for poinng this out. We fully agree that the accuracy (such as root-mean-square error of energy and forces) of the ML potenal does not necessarily guarantee the stability of larger length/me scales MD simulaons. In addion to the acve learning process, we monitored the extrapolaon grade value during the long-mescale MD simulaons of our large system. Specifically, configuraons with an extrapolaon grade value γ ≤ 1 are considered interpolaon structures, whereas those with γ > 1 are regarded as extrapolaon structures. Because the extrapolaon grade value is posively correlated with the errors in forces and energies,30 it allows us to effecvely monitor and ensure the accuracy of our MD simulaons. The extrapolaon grade values obtained during the 2 ns NVE ensemble MD simulaons are provided in SI (secon Extrapolaon grade of SI, Fig. S11-S29).\nIn accordance with the reviewers’ recommendations, we extracted structural snapshots every 5 fs from the 2 ns NVE MD simulation trajectory, resulting in a total of 400,000 configurations. These configurations were used to analyze the oxygen and hydrogen atomic probability densities (Fig. S49–S67) and the water density profiles along the Z axis (Fig. S68–S86). Additionally, snapshots were extracted every 100 fs to obtain 20,000 configurations for calculating the O–O and O–H radial distribution functions (Fig. S87–S90). The resulting oxygen and hydrogen atomic probability densities are consistent with those previously reported for suspended monolayer graphene systems by Rashmi et al.25 and Wang et al..31 The corresponding clarificaon has been added on page 8 of the manuscript and is highlighted with a background color.\n1. Are the MB-mu and MB-alpha, which may be applicable for the water/air interface, applicable for water/solid interfaces? I mean, do the dipole moment surface and polarizability surface cover the water structure modulated by the solid surface? Such a validaon would be required.\nWe appreciate the reviewer’s comment. In our work, the ACE potenal was trained using energies and forces calculated from density funconal theory, which fully accounts for the effects of both graphene and the CaF₂ substrate. During the MD simulaons, we employed the ACE potenal, which inherently includes the influence of graphene and the substrate’s polarizable effects on the orientaon of interfacial water molecules. Prior to the vSFG spectra simulaons, the dipole moments and polarizabilies of water molecules were calculated using the MB-μ and MB-α potenals, respecvely., which were parameterized on a broad set of gas-phase and small cluster geometries.22 We acknowledge that the MB-μ and MB-α potenals do not explicitly incorporate the induced dipole moments or atomic polarizabilies in water molecules that arise from interacons with graphene or substrate atoms. Instead, the dipole moment and polarizability of each water molecule depend solely on its geometry and orientaon obtained from the trajectories generated by our ACE potenal. This strategy—considering only the geometry and orientaon of water molecules for vSFG spectra simulaons—has been widely adopted by other research groups studying solid/liquid interfaces.11, 31-35 Importantly, the vSFG spectra obtained using this approach show good agreement with both experimental measurements and previous computaonal results for the suspended graphene/water11, 31, 32, 36 and CaF₂/water1, 11, 16, 37 interfaces. Furthermore, Rashmi et al.25 recently performed the path-integral molecular dynamics simulaons of the suspended graphene/water system and computed its vSFG spectra using the MB-μ and MB-α potenals, demonstrang good agreement with experimental data and yielding results consistent with our simulaons. This consistency further supports that our current strategy provides a reasonable and reliable descripon of the interfacial systems invesgated in this study. To clarify this point, 24 we have added the following sentences on pages 26-27 of the manuscript, highlighted with a colored background.\n“We acknowledge that the MB-µ and MB-α do not explicitly incorporate the induced dipole moments or atomic polarizabilies in water molecules that arise from interacons with graphene or substrate atoms. Instead, the dipole moment and polarizability of each water molecule depend solely on the geometry and orientaon of water molecules obtained from the trajectories generated by our ACE potenal. This strategy—considering only the geometry and orientaon of water molecules for vSFG spectra simulaons—has been widely adopted by other research groups studying solid/liquid interfaces and shows good agreement with the experimental measurements.[55, 57, 68, 86, 153–155]”\n2. The SFG spectrum for the water/air interface looks a bit deviated from a spectrum computed similarly (J. Phys. Chem. B 122, 4356 (2018).). The presented spectrum (Fig. 2c), in parcular, the feature of the 3600 shoulder peak, looks like the low-temperature spectrum (238K) in the reference. Is it because the MLIP is learned from water/graphene interface?\nWe appreciate the reviewer’s comment. It is true that in the referenced study (J. Phys. Chem. B 122, 4356 (2018)), the shoulder near 3600 cm⁻¹ in the simulated spectra appears at low temperatures and gradually weakens as the temperature increases, disappearing completely at room temperature and above. However, several experimentally measured heterodyne-detected vSFG spectra, as shown in Figure 4a–d, obtained by different research groups, clearly exhibit a shoulder feature near this frequency at or close to room temperature for both H₂O and D₂O at the water/air interface. In the study corresponding to Figure 4b in Ref. 39, the authors atributed the broad shoulder (indicated by the black arrow) to intramolecular vibraonal coupling between OD chromophores. This feature is suppressed in the uncoupled case of dilute HOD in H₂O.\nIn the study corresponding to FIG. 1 in Ref 28, the authors compared vSFG spectra calculated using the surface-specific velocity–velocity correlaon funcon (ssVVCF) and 25 the dipole moment-polarizability (μ–α) correlaon funcon. They emphasized that the shoulder appears only in the μ–α correlaon funcon spectrum. This occurs because the ssVVCF approach projects the transion dipole and polarizability moments onto atomic velocies, which cannot capture key spectral features arising from vibraonal coupling—hence the absence of the shoulder in ssVVCF spectra. Furthermore, Litman et al.38 computed vSFG spectra using molecular dynamics trajectories generated with various funconals (PBE-D3, BLYP-D3, revPBE-D3, HSE06, B3LYP, and revPBE0-D3). In all cases, the simulated spectra exhibited a shoulder near 3600 cm⁻¹ (Figure 2 in Ref. 38).\nTaken together, both experimental observaons and μ–α correlaon funcon–based simulaons consistently display a shoulder feature at or near room temperature. Therefore, we believe that our simulated vSFG spectra for the water/air interface more accurately reproduce the experimentally observed spectral characteriscs. An explanaon regarding the origin of the shoulder feature in the vSFG spectrum of the water/air interface has been included on page 10 of the manuscript.\nFigure 4 a, b, c, and d, Heterodyne-detected vSFG spectra (Imχ(2)) from Ref.40(CC BY 3.0), Ref.31(CC BY), Ref.41(CC BY 4.0) and Ref.42(CC BY-NC 4.0), respecvely.\n26\n3. If the authors would like to suggest that the experimentally observed vSFG signal arose from signal cancellaon, it would be beter to compare computed signals with some measured signals and discuss.\nWe appreciate the reviewer’s comment. As shown in Figure S13 1L in Ref. 1 and Figure 5a, independently measured vSFG spectra from two research groups for monolayer graphene on a CaF₂ substrate do not exhibit a disnct peak near 3600 cm⁻¹.1, 37 Although a weak shoulder feature can be observed around 3600 cm⁻¹, as indicated by the blue arrow in Figure 5a, this suggests that the monolayer graphene supported on CaF₂ exhibits a hydrophilic character. In contrast, for four-layer graphene on a CaF₂ substrate (Figure S13 4L in Ref. 1), an obvious peak (or shoulder) appears around 3600 cm⁻¹, which corresponds to dangling OH groups poinng toward the graphene surface. This observaon indicates that the four-layer graphene supported on CaF₂ behaves more hydrophobically.\nFigure 5 a, Raw vSFG spectrum of monolayer graphene on a CaF₂ substrate, reproduced from Ref.37(CC BY). b, The squared modulus |χ(2)|² of the vSFG spectra (sum of the squared imaginary and real components) for the S L and S L 64 1 64 4 systems obtained from our simulaons.\nFigure 5b presents our simulated vSFG spectra obtained with the same number of intercalated water molecules. For the monolayer graphene/CaF₂ system, with some intercalated water, the peak near 3600 cm⁻¹ remains weak, similar to the experimental spectra shown in Figure S13 1L in Ref. 1 and Figure 5a. With some intercalated water, the simulated spectrum of the four-layer graphene/CaF₂ system shows only a weak feature 27 near 3600 cm⁻¹, whereas the experimental measurements display a more pronounced peak (or shoulder) at this frequency.\nThis discrepancy movated us to calculate the free energy of intercalated water molecules relave to those water molecules in the slab region. Our free-energy results reveal that water intercalaon is thermodynamically favorable for monolayer graphene/CaF₂ but unfavorable for four-layer graphene/CaF₂. Even when the same number of intercalated water molecules is introduced in both the monolayer graphene/CaF₂ and four-layer graphene/CaF₂ systems, the simulated vSFG spectra show similar peak intensies around 3600 cm⁻¹. However, since water intercalaon is thermodynamically unfavorable for four-layer graphene/CaF₂, the presence of intercalated water molecules in this system is unlikely. This finding explains the presence of a clear 3600 cm⁻¹ peak (shoulder) in the experimental vSFG spectrum of the four-layer graphene/CaF₂ system but not in that of the monolayer graphene/CaF₂ system. To clarify this point, we have added the following sentences on page 18 of the manuscript, highlighted with a colored background.\n“Independently measured vSFG spectra from different research groups for monolayer graphene on a CaF substrate do not exhibit a disnct peak near 3665 2 cm−1.[35, 54, 58] Although a weak shoulder feature around 3650 cm−1 was observed in the vSFG spectra measured by Dreier et al.,[58] this suggests that the monolayer graphene supported on CaF exhibits a hydrophilic character. In ` contrast, for graphene films thicker than three layers, a disnct peak (or shoulder) emerges in the frequency range of 3600-3700 cm−1,[35] whereas this feature is absent in systems with one to three graphene layers.”\n4. The mestep of the DFT-MD simulaon (1.0 fs) is a bit aggressive for systems containing O-H stretch modes. It is beter to show the total energies of DFT-MD in addion to those of MLIP-MD. Also, how long DFT-MD simulaons are performed is not clearly menoned.\nWe appreciate the reviewer’s comment. We agree that the 1.0 fs mestep is close to the upper limit for systems containing O–H stretching vibraons. In our workflow, the 1.0 fs mestep was used only during the short (tens to hundreds picoseconds, Table S1) on-the- 28 fly machine-learning force field molecular dynamics (OFMLFFMD) simulaons to generate part of the training dataset. These trajectories were not used for interfacial structure analysis or vSFG spectral simulaon. During the acve learning molecular dynamics (ALMD) simulaon, which idenfies more representave structures on the potenal energy surface for refining the potenal, a smaller mestep of 0.5 fs was employed. All equilibraon and producon MD simulaons for large systems used for interfacial structure analysis and vSFG spectral calculaons were performed with the final ACE potenal using a mestep of 0.5 fs. Detailed informaon about the OFMLFFMD and ALMD simulaons is provided in Table S1 of the SI.\nThe following sentences have been added on page 24 of the manuscript:\n“Although a 1.0 fs mestep, close to the upper limit for systems containing O–H stretching vibraons, was used in these short MD simulaons, the resulng trajectories were employed only for training data generaon and were not ulized for interfacial structure analysis or vSFG spectral simulaons.”\nAs described in Secon 5.1, we did not perform AIMD simulaons. Instead, we used the on-the-fly machine-learning force field implemented in VASP to carry out the MD simulaons. This strategy markedly reduces the computaonal cost and thereby enables efficient sampling of the configuraonal space for our complex system. To validate the accuracy of our ACE potenal, we performed MD simulaons with LAMMPS using our ACE potenal for several systems with cell sizes that could directly do DFT energy calculaons. Each simulaon lasted 2 ns, including 1 ns of equilibraon in the NVT ensemble using the canonical sampling through velocity rescaling (CSVR) thermostat, followed by 1 ns in the NVE ensemble. The mestep for these simulaons was 0.5 fs. The model used for validaon is shown in Figure 6. From the 1 ns NVE ensemble MD trajectories, one structure was extracted every 5 ps, resulng in a total of 200 configuraons. Single-point DFT (BLYP-D4) energy calculaons were subsequently performed for these configuraons. The potenal energies predicted by the ACE potenal were then compared to the corresponding DFT energies, as shown in Figure 6.\n5. I think that the effect of the cross-correlation term on the spectra is important and should be presented in the supplementary materials.\nWe appreciate the reviewer’s comment. Since the length of our trajectory (2 ns) is sufficiently long, the cross-correlaon term, which accounts for correlaons between the dipole moment and polarizability of neighboring molecules within a defined cutoff distance, was explicitly included in our vSFG spectral simulaons. Specifically, only 31 contribuons from water molecule pairs with an oxygen–oxygen distance shorter than 4.0 Å were considered, as described in Equaon (2) of the manuscript. This cutoff value has been adopted in several previous studies.23-25, 43\nTo clarify the effect of the cross-correlaon term on the vSFG spectra, we selected four systems: suspended monolayer graphene (L ), substrate-supported monolayer graphene 1 (SL ), and structures with intercalated water between the substrate and graphene layers 1 of varying thickness (S L and S L ). For the Im(χ(2)) spectra of L and SL (Figure 7a), the 64 1 64 4 1 1 cross-correlaon term has a negligible effect on the high-frequency dangling OH peak but significantly affects the broad low-frequency peak, resulng in a reduced peak amplitude and a blue shi of approximately 100 cm⁻¹. This effect is similarly observed in the squared modulus |χ(2)|²(Figure 7b), where neglecng the cross-correlaon term leads to a decreased intensity of the broad low-frequency peak and a blue shi, but a relavely higher intensity of the high-frequency dangling OH peak.\nFor the water-intercalated systems (S L and S L ), omission of the cross-correlaon term 64 1 64 4 results in lower intensies for both the high-frequency dangling OH peak and the low-frequency hydrogen-bonded OH peak in the Im(χ(2)) spectra (Figure 7c), with a more pronounced effect on the low-frequency peak. Consistent with the non-intercalated systems, excluding the cross-correlaon term causes a blue shi of about 100 cm⁻¹ in the low-frequency hydrogen-bonded OH peak. The squared modulus |χ(2)|² for S L and S L 64 1 64 4 (Figure 7d) follows a similar trend to that of L and SL : neglecng the cross-correlaon 1 1 term decreases the intensity of the broad low-frequency peak with a blue shi, while relavely increasing the intensity of the high-frequency dangling OH peak.\nTo understand why the cross-correlaon term affects the high-frequency peak of the intercalated systems (S L and S L ) more strongly than in L and SL , we calculated the 64 1 64 4 1 1 vSFG spectra of the intercalated water molecules alone, with and without the cross-correlaon term. As shown in Figure 7e, removing the cross-correlaon term enhances the negave high-frequency peak in both S L1 and S L , an effect that is even more 64 64 4 pronounced in the squared modulus |χ(2)|² (Figure 7f). These results demonstrate that the 32 cross-correlaon contribuon is essenal for accurately describing the vSFG response of water confined in restricted spaces. The above clarificaon and Figure 7, which discuss the effect of the cross-correlaon term in the me-correlaon funcon, have been added to pages 93–94 of the SI.\nFigure 7 vSFG spectra calculated with and without the cross-correlaon term in the me-correlaon funcon. Systems include floang monolayer graphene (L ), substrate-supported monolayer graphene (SL ), and structures with 1 1 intercalated water between substrate and graphene layers of varying thickness (S 64 L 1 and S 64 L 4 ). Figures a and b correspond to L 1 and SL 1 , whereas figures c and d correspond to S 64 L 1 and S 64 L 4 . Figures a and c show the imaginary part, Im (χ(2)), and figures b and d show the squared modulus, |χ(2)|2. Figures e and f present Im (χ(2)) and |χ(2)|2, respecvely, 33 for the intercalated water molecules in S L and S L . Solid lines represent spectra including both auto- and cross- 64 1 64 4 correlaon terms; dashed lines include only the autocorrelaon term.\n6. I noced that the citaon (ref. 72) is not consistent.\nWe appreciate your careful observaon. We have corrected the inconsistency in the citaon of this reference, which is now listed as reference 70 (Front. Phys. 9:740627). In addion, we have carefully rechecked all other references to ensure their accuracy and consistency throughout the manuscript.\n7. Some labels overlap in Fig. 6a.\nThank you for poinng this out. In response, we have adjusted Figure 6a to remove the overlapping labels and improve the overall readability of the figure. The updated figure is shown below:\nFigure 8 a, Free energy of intercalated water in monolayer and four-layer graphene systems at different coverages. The free energy per water molecule in the corresponding water slab is set as the reference (zero). The plot shows the relave free energy of each intercalated water molecule compared to the reference. b, Schemac illustraon showing that water molecules energecally prefer intercalaon between monolayer graphene and substrate over surface adsorpon. c, Schemac illustraon showing that water molecules energecally prefer surface adsorpon on four-layer graphene over intercalaon.\n(1) Kim, D.; Kim, E.; Park, S.; Kim, S.; Min, B. K.; Yoon, H. J.; Kwak, K.; Cho, M. Wetability of graphene and interfacial water structure. Chem 2021, 7 (6), 1602-1614.\n(2) Berry, V. Impermeability of graphene and its applicaons. Carbon 2013, 62, 1-10.\n(3) Arif, T.; Colas, G.; Filleter, T. Effect of Humidity and Water Intercalaon on the Tribological Behavior of Graphene and Graphene Oxide. ACS Appl. Mater. Interfaces 2018, 10 (26), 22537-22544.\n(4) Xu, K.; Cao, P.; Heath, J. R. Graphene Visualizes the First Water Adlayers on Mica at Ambient Condions. Science 2010, 329 (5996), 1188-1191.\n(5) Severin, N.; Lange, P.; Sokolov, I. M.; Rabe, J. P. Reversible Deweng of a Molecularly Thin Fluid Water Film in a So Graphene–Mica Slit Pore. Nano Lett. 2012, 12 (2), 774-779.\n(6) Rezania, B.; Dorn, M.; Severin, N.; Rabe, J. P. Influence of graphene exfoliaon on the properes of water-containing adlayers visualized by graphenes and scanning force microscopy. J. Colloid Interface Sci. 2013, 407, 500-504.\n(7) Lee, M. J.; Choi, J. S.; Kim, J.-S.; Byun, I.-S.; Lee, D. H.; Ryu, S.; Lee, C.; Park, B. H. Characteriscs and effects of diffused water between graphene and a SiO2 substrate. Nano Res. 2012, 5 (10), 710-717.\n(8) Lee, D.; Ahn, G.; Ryu, S. Two-Dimensional Water Diffusion at a Graphene–Silica Interface. J. Am. Chem. Soc. 2014, 136 (18), 6634-6642.\n(9) Temmen, M.; Ochedowski, O.; Schleberger, M.; Reichling, M.; Bollmann, T. R. J. Hydraon layers trapped between graphene and a hydrophilic substrate. New J. Phys. 2014, 16 (5), 053039.\n(10) Kim, J.-S.; Choi, J. S.; Lee, M. J.; Park, B. H.; Bukhvalov, D.; Son, Y.-W.; Yoon, D.; Cheong, H.; Yun, J.-N.; Jung, Y.; et al. Between Scylla and Charybdis: Hydrophobic Graphene-Guided Water Diffusion on Hydrophilic Substrates. Sci. Rep. 2013, 3 (1), 2309.\n(11) Wang, Y.; Tang, F.; Yu, X.; Chiang, K.-Y.; Yu, C.-C.; Ohto, T.; Chen, Y.; Nagata, Y.; Bonn, M. Interfaces govern the structure of angstrom-scale confined water soluons. Nat. Commun. 2025, 16 (1), 7288.\n(12) Cen, J.; Liang, F.; Chen, D.-L.; Zhang, L.; Yang, N.; Zhu, W. Adsorpon of Water Molecule on Calcium Fluoride and Magnesium Fluoride Surfaces: A Combined Theorecal and Experimental Study. J. Phys. Chem. C 2020, 124 (14), 7853-7859.\n(13) Koma, A.; Saiki, K.; Sato, Y. Heteroepitaxy of a two-dimensional material on a three-dimensional material. Appl. Surf. Sci. 1990, 41-42, 451-456.\n(14) Song, H.; Zhou, F.; Yan, S.; Su, X.; Wu, H.; Wu, Q.; Gao, Y.; Chen, R.; Chen, T.; Yao, J.; et al. Enhanced Transport and Optoelectronic Properes of van der Waals Materials on CaF2 Films. Nano Lett. 2023, 23 (11), 4983-4990.\n(15) Wang, Y.; Seki, T.; Liu, X.; Yu, X.; Yu, C.-C.; Domke, K. F.; Hunger, J.; Koper, M. T. M.; Chen, Y.; Nagata, Y.; et al. Direct Probe of Electrochemical Pseudocapacive pH Jump at a Graphene Electrode**. Angew. Chem., Int. Ed. 2023, 62 (10), e202216604.\n(16) Lesnicki, D.; Zhang, Z.; Bonn, M.; Sulpizi, M.; Backus, E. H. G. Surface Charges at the CaF2/Water Interface Allow Very Fast Intermolecular Vibraonal-Energy Transfer. Angew. Chem., Int. Ed. 2020, 59 (31), 13116-13121.\n(17) Khab, R.; Backus, E. H. G.; Bonn, M.; Perez-Haro, M.-J.; Gaigeot, M.-P.; Sulpizi, M. Water orientaon and hydrogen-bond structure at the fluorite/water interface. Sci. Rep. 2016, 6 (1), 24287.\n(18) Lysogorskiy, Y.; Oord, C. v. d.; Bochkarev, A.; Menon, S.; Rinaldi, M.; Hammerschmidt, Hammerschmidt, T.; Mrovec M Thompson A Csányi G Ortner C et al Performant implementaon of the atomic cluster expansion (PACE) and applicaon to copper and silicon. npj Comput. Mater. 2021, 7 (1), 97.\n(19) Bochkarev, A.; Lysogorskiy, Y.; Menon, S.; Qamar, M.; Mrovec, M.; Drautz, R. Efficient parametrizaon of the atomic cluster expansion. Phys. Rev. Mater. 2022, 6 (1), 013804.\n(20) Drautz, R. Atomic cluster expansion for accurate and transferable interatomic potenals. Phys. Rev. B 2019, 99 (1), 014104.\n(21) Medders, G. R.; Paesani, F. Many-Body Convergence of the Electrostac Properes of Water. J. Chem. Theory Comput. 2013, 9 (11), 4844-4852.\n(22) Medders, G. R.; Paesani, F. Infrared and Raman Spectroscopy of Liquid Water through “First-Principles” Many-Body Molecular Dynamics. J. Chem. Theory Comput. 2015, 11 (3), 1145-1154.\n(23) Medders, G. R.; Paesani, F. Dissecng the Molecular Structure of the Air/Water Interface from Quantum Simulaons of the Sum-Frequency Generaon Spectrum. J. Am. Chem. Soc. 2016, 138 (11), 3912-3919.\n(24) Moberg, D. R.; Straight, S. C.; Paesani, F. Temperature Dependence of the Air/Water Interface Revealed by Polarizaon Sensive Sum-Frequency Generaon Spectroscopy. J. Phys. Chem. B 2018, 122, 4356.\n(25) Rashmi, R.; Balogun, T. O.; Azom, G.; Agnew, H.; Kumar, R.; Paesani, F. Revealing the Water Structure at Neutral and Charged Graphene/Water Interfaces through Quantum Simulaons of Sum Frequency Generaon Spectra. ACS Nano 2025, 19 (4), 4876-4886.\n(26) Wang, Y.; Seki, T.; Yu, X.; Yu, C.-C.; Chiang, K.-Y.; Domke, K. F.; Hunger, J.; Chen, Y.; Nagata, Y.; Bonn, M. Chemistry governs water organizaon at a graphene electrode. Nature 2023, 615 (7950), E1-E2.\n(27) Brandenburg, J. G.; Zen, A.; Alfè, D.; Michaelides, A. Interacon between water and carbon nanostructures: How good are current density funconal approximaons? J. Chem. Phys. 2019, 151 (16), 164702.\n(28) Ohto, T.; Usui, K.; Hasegawa, T.; Bonn, M.; Nagata, Y. Toward ab inio molecular dynamics modeling for sum-frequency generaon spectra; an efficient algorithm based on surface-specific velocity-velocity correlaon funcon. J. Chem. Phys. 2015, 143 (12), 124702.\n(29) Tang, F.; Ohto, T.; Sun, S.; Rouxel, J. R.; Imoto, S.; Backus, E. H. G.; Mukamel, S.; Bonn, M.; Nagata, Y. Molecular Structure and Modeling of Water–Air and Ice–Air Interfaces Monitored by Sum-Frequency Generaon. Chem. Rev. 2020, 120 (8), 3633-3667.\n(30) Lysogorskiy, Y.; Bochkarev, A.; Mrovec, M.; Drautz, R. Acve learning strategies for atomic cluster expansion models. Phys. Rev. Mater. 2023, 7 (4), 043801.\n(31) Wang, Y.; Tang, F.; Yu, X.; Ohto, T.; Nagata, Y.; Bonn, M. Heterodyne-Detected Sum-Frequency Generaon Vibraonal Spectroscopy Reveals Aqueous Molecular Structure at the Suspended Graphene/Water Interface. Angew. Chem., Int. Ed. 2024, 63 (20), e202319503.\n(32) Ohto, T.; Tada, H.; Nagata, Y. Structure and dynamics of water at water–graphene and water–hexagonal boron-nitride sheet interfaces revealed by ab inio sum-frequency generaon spectroscopy. Phys. Chem. Chem. Phys. 2018, 20 (18), 12979-12985.\n(33) Azom, G.; Milet, A.; David, R.; Kumar, R. From Graphene Oxide to Graphene: Changes in Interfacial Water Structure and Reacvity Using Deep Neural Network Force Fields. J. Phys. Chem. C 2024, 128 (39), 16437-16453.\n(34) Tabassum, T.; Das, B.; Park, C.; Ghosh, S.; Naujoks, D.; Piontek, S. M.; Ludwig, A.; Marx, D.; Petersen, P. B. Hydrophobic Versus Hydrophilic Nature of the Gold-Water Interface Determined by Fluctuang Local Water Orientaon. Phys. Rev. Lett. 2025, 135 (11), 116201.\n(35) Kroul, O.; Pezzo, S.; Gaigeot, M.-P.; Předota, M. Phase-Sensive Vibraonal SFG Spectra from Simple Classical Force Field Molecular Dynamics Simulaons. J. Phys. Chem. C 2020, 124 (28), 15253-15263.\n(36) Xu, Y.; Ma, Y.-B.; Gu, F.; Yang, S.-S.; Tian, C.-S. Structure evoluon at the gate-tunable suspended graphene–water interface. Nature 2023, 621 (7979), 506-510.\n(37) Dreier, L. B.; Liu, Z.; Narita, A.; van Zadel, M. J.; Müllen, K.; Tielrooij, K. J.; Backus, E. H. G.; Bonn, M. Surface-Specific Spectroscopy of Water at a Potenostacally Controlled Supported Graphene Monolayer. J. Phys. Chem. C 2019, 123 (39), 24031-24038.\n(38) Litman, Y.; Lan, J.; Nagata, Y.; Wilkins, D. M. Fully First-Principles Surface Spectroscopy with Machine Learning. J. Phys. Chem. Lett. 2023, 14 (36), 8175-8182.\n(39) Sopkin, I. V.; Weeraman, C.; Pieniazek, P. A.; Shalhout, F. Y.; Skinner, J. L.; Benderskii, A. V. Hydrogen bonding at the water surface revealed by isotopic diluon spectroscopy. Nature 2011, 474 (7350), 192-195.\n(40) Nagata, Y.; Hasegawa, T.; Backus, E. H. G.; Usui, K.; Yoshimune, S.; Ohto, T.; Bonn, M. The surface roughness, but not the water molecular orientaon varies with temperature at the water–air interface. Phys. Chem. Chem. Phys. 2015, 17 (36), 23559-23564, 10.1039/C5CP04022A.\n(41) Seki, T.; Yu, C.-C.; Chiang, K.-Y.; Greco, A.; Yu, X.; Matsumura, F.; Bonn, M.; Nagata, Y. Ions Speciaon at the Water–Air Interface. J. Am. Chem. Soc. 2023, 145 (19), 10622-10630.\n(42) Chiang, K.-Y.; Hunger, J.; Bonn, M.; Nagata, Y. Experimental quanficaon of nuclear quantum effects on the hydrogen bond of liquid water. Sci. Adv. 11 (14), eadv7218.\n(43) Rashmi, R.; Paesani, F. Dissecng the Molecular Structure of the Air/Ice Interface from Quantum Simulaons of the Sum-Frequency Generaon Spectrum. J. Am. Chem. Soc. 2025, 147 (2), 1903-1910."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I have reviewed the authors' comprehensive point-by-point responses to my comments and the corresponding revisions made to the manuscript and supplementary information.\nThe authors have convincingly justified their choice of employing ACE potentials. The explicit discussion added to Section 2.2 and the provision of all training datasets, potentials, and simulation files on Zenodo greatly enhance the transparency, reproducibility, and overall robustness of the study. The assumptions regarding the absence of interlayer water intercalation in multilayer graphene have been properly contextualized with experimental evidence and clearly stated as model boundaries. The citation error has been corrected, and the discussion on kinetic versus thermodynamic control of intercalation in multilayer graphene is a valuable addition. The rationale for selecting the CaF2(111) surface is now clearly explained, and the expected limitations of the trained potentials concerning water coverage are appropriately documented. The brief explanation of the vSFG calculation methodology is also helpful for readers.\nOverall, the authors have thoughtfully addressed all my initial comments. The manuscript has been substantially improved in clarity, depth, and scholarly rigor. I recommend the manuscript for publication in its present form.", "reviewer_2": "The authors have made a clear effort to address the reviewers’ concerns and to improve both the positioning and technical transparency of the work. The central result that pristine graphene is intrinsically hydrophobic, and that “hydrophilic” vSFG signatures reported for substrate-supported monolayer graphene can arise from the coexistence of adsorbed and intercalated water (leading to partial cancellation of spectral features) is now presented more clearly and connected more explicitly to the experimental context.\nThere are two remaining points that may require some minor clarifications:\n1) The manuscript now explains why NQEs were not included and cites prior work indicating that NQEs mainly induce a red shift in the OH-stretch region without large changes in the overall spectral profile. Since the present paper also places weight on intercalation thermodynamics (free energies) in addition to spectral trends, it would be useful to add one brief cautionary sentence noting that NQEs could, in principle, affect free-energy differences for confined/intercalated water even when the spectral line shape is only modestly altered. This can be framed as a current limitation of the model and a direction for future work.\n2) The rationale for choosing BLYP-D4 is now clearer. Benchmarking studies of water–carbon interactions have identified BLYP-D4 as one of the best-performing GGA-level approaches for this class of systems (within the set of tested GGAs). However, it is also well documented that BLYP can provide an imperfect description of bulk liquid water in ab initio simulations (e.g., yielding an overly structured liquid and strongly underestimated diffusivity compared to experiment). More generally, water properties can be sensitive to the chosen electronic-structure reference: for example, Galib et al. [J. Chem. Phys. 146, 244501 (2017)] explicitly compared results obtained with dispersion-corrected revPBE and BLYP, highlighting the dependence of predicted observables on the underlying interaction model. More recently, de Hijes et al. [J. Chem. Phys. 161, 131102 (2024)] reported extensive comparisons of ML potentials of water derived from various DFT approximations and concluded that further improvements in DFT approximations are necessary to enable accurate quantitative predictions for water.\nThese studies, together with the extensive literature on DFT-based simulations of water, suggest that, beyond the error compensation discussed by the authors in connection with neglected NQEs, additional cancellation of errors between graphene–water and water–water interactions may occur at the BLYP-D4 level. This possibility implies that the interpretation of the vSFG response may be less robust than one would ideally wish. It would therefore be important to state explicitly that the conclusions are drawn at the BLYP-D4 reference level and that functional sensitivity remains an inherent limitation of DFT-based simulations of water.\nWith the above clarifications, the revised manuscript is technically sound and significantly improved in presentation and transparency.", "reviewer_3": "The authors conducted additional analysis and simulations required by the reviewers, which improved the quality of the manuscript significantly. The revised manuscript contains insights into interfacial water structure at graphene on the basis of sound molecular dynamics simulations. Now I recommend to accept this manuscript for publication."}, "author_response": ""}]

Decision letter