Reflection on the horizontal and vertical organization of water at oxide interfaces and their SFG vibrational signatures

Abstract

The hydrophilic and hydrophobic wetting behaviors of SiO₂ and Al₂O₃ surfaces at PZC conditions are examined at the molecular level using graph theory. Hydrophobic wettability is manifested by a remarkable, collectively organized horizontal 2D hydrogen-bonded network (2DN) of water, anchored to the oxide surface through a dense array of vertically oriented hydrogen-bonded motifs. In contrast, hydrophilic wettability is characterized solely by vertical hydrogen-bonded motifs, similar to those that anchor the 2DN at hydrophobic aqueous interfaces. Analysis of these vertical motifs – considering both their structural features and SFG spectroscopic signatures – provides a direct molecular-level connection between the microscopic organization of interfacial water and the macroscopic manifestations of hydrophilic and hydrophobic behavior.

---

1 Introduction

The structure that water can adopt when it is in contact with a solid surface is of paramount importance for the chemical and physical processes that can occur at aqueous interfaces, in, e.g., geochemistry, atmospheric science, electrochemistry, and corrosion, to cite a few fields.¹ Typically, heterogeneous catalysis, electrochemistry and energy storage are central fields in chemistry where the structure of water at the interface with a solid plays a pivotal role in the outcome of chemical reactions and charge transfers within the electric double layer. Transport of pollutants in underground water is another area in research and technology in which the structure of interfacial water at the interface with minerals plays a crucial role. Water transport in nanofluidics is also controlled by the structure and interactions of water at the interface with solids.

A lot of research is being conducted to figure out the structure of aqueous interfaces, both experimentally and theoretically.^1–14 Macroscopically, the wetting at aqueous interfaces is classified as hydrophilic or hydrophobic, often based on contact angle measurements.^2,15 From the point of view of water, wetting at hydrophobic surfaces can be defined as water avoiding as much as possible the contact with the surface, which, at the molecular level, can be translated into an inability or lack of ability to form hydrogen bonds between the surface and the water. In contrast, hydrophilic wetting can be macroscopically defined as water willing to interact with the surface, which, at the molecular level, is nothing more than the high ability of water to form hydrogen bonds with the surface.

From a structural point of view, hydrogen bonding is hence central to the molecular concepts of hydrophilic and hydrophobic wetting.^1,2 While the knowledge of hydrogen bonding can be directly obtained from molecular dynamics simulations of aqueous interfaces, it is not as straightforward to extract such knowledge, directly, from microscopic experimental measurements at aqueous interfaces. Surface sensitive Sum Frequency Generation (SFG) spectroscopy^16–18 is certainly the technique of choice for obtaining information on water wetting surfaces, at the molecular level. In this work, our interest lies on aqueous oxides for which numerous SFG experiments have been conducted,^1,17–24 in particular on the prototype silica oxide in contact with liquid water in various pH and electrolytic concentration conditions.¹ Heterodyne HD-SFG experiments have provided direct information on the orientation of water in contact with the electric double layer of the oxide surface (from the knowledge of the positive/negative sign of the spectroscopic bands) and on the strength of the ‘hydrophilic’ surface–water interactions (from the knowledge of the positions of the bands). However, this information on the structure of water at the interface with the oxide surface is still fragmentary, as it does not yet allow to understand how water wets the oxide, i.e., what the actual collective microscopic structure of water at the interface with the surface is. Only this information will give us a detailed understanding on the hydrophilic vs. hydrophobic wetting of the surface oxides.

At the macroscopic level, oxides are usually classified as hydrophilic. At the molecular level, this ability of water to form hydrogen bonds with the oxide surface has been observed on several occasions through the positions of the SFG bands in the range ∼3400–3000 cm⁻¹, which corresponds to moderate to strong hydrogen bonding forces.^1,25 On the other hand, the hydrophobic character of the water at the interface with an oxide has also been measured and documented in SFG. This is achieved through the appearance of the spectroscopic signature of free OH groups of water around 3600–3800 cm⁻¹.^26–29 This has especially been documented by some of us in ref. 29 in the case of silica oxides. At aqueous oxides, this free OH spectroscopic feature is generally accompanied by the spectroscopic signatures of hydrogen bonded OH, see for example ref. 29, which makes it difficult to conclude on the global hydrophobic–hydrophilic nature of interfacial water. One generally speaks in terms of hydrophobic areas or patches over the oxide surface that induce the presence of hydrophobic water patches, see e.g. ref. 29 in the case of silica oxide. Unless the SFG signal solely comprises the free OH band, there is still no definite quantification of hydrophobicity of the aqueous oxide interface obtained from SFG measurements.

Several microscopic descriptors of hydrophobicity have been developed in the last decades through MD simulations, based on density fluctuations,^30–32 bond interactions,³³ or orientation probes³⁴ of interfacial water. With the hydrogen bonding ability of water in mind, we developed in ref. 15 a molecular descriptor based on the balance of surface–water and water–water H-bonds, denoted x_H: it measures the ratio of horizontal vs. vertical H-bonds that the water can form at an interface. Horizontal H-bonds are the ones made in between the water molecules that are located within the interfacial layer (the Binding Interfacial Layer, BIL, as defined in ref. 35–37), while vertical H-bonds are the ones formed between water in the BIL and the top-surface sites. Section 2.2 of this paper will describe the concepts in more detail. These two H-bond types can be easily extracted and counted from molecular dynamics simulations of aqueous interfaces. The molecular x_H descriptor has been shown in ref. 15 to correlate with macroscopic measurements of contact angles and has been shown to rationalize SFG macroscopic measurements into hydrophilicity and hydrophobicity wettability characters of water on model interfacial systems.

The balance between horizontal and vertical H-bonds at aqueous interfaces has been further correlated with solvation free energies,^38,39 hence correlating this descriptor with the ability of aqueous interfaces to drive adsorption/desorption of solutes as well as their ability to adapt to surface fields. The same authors have also applied these concepts to water at the interface with biomolecules,⁴⁰ showing that the balance between horizontal and vertical H-bonds correlates linearly with the intensity of THz spectroscopic features of water H-bond stretching and librational inter-molecular motions.

In this work, the molecular descriptor x_H is associated with algorithmic graph theory and theoretical SFG spectroscopy in order to analyze the ability of water to form horizontal and vertical hydrogen bonds at the interface with silica and alumina oxides, and to provide a direct measurement of the hydrophilic and hydrophobic characteristics of these aqueous interfaces. SFG signatures are especially dissected in order to possibly find relationships, at the molecular level, between the spectroscopic signatures and the hydrophilic–hydrophobic properties of these aqueous interfaces. Along this theoretical journey, we introduce surface–water H-bonded vertical motifs, directly obtained from graph theory, and their individual spectroscopic signatures that could ultimately constitute a database of motif structure – the SFG signature. The analyses of the molecular dynamics trajectories are conducted through the concepts of graph theory,^41–43 that will in particular be shown to be essential for the characterization of the collective hydrophobic organization of water at silica and alumina aqueous oxides.

Our paper is hence organized as follows. Section 2 presents the whole methodology applied throughout the work: DFT-MD simulations, the x_H molecular descriptor of hydrophilicity/hydrophobicity, the description of the molecular graphs, including the graph algorithms for their analysis in terms of collective components (for the horizontal organization of water) and of sub-graphs (for the vertical organization of water), and finishing with SFG spectroscopic calculations.

Section 3 provides the results of the silica and alumina oxide facets that behave hydrophobic or hydrophilic at the interface with liquid water, from the molecular descriptor x_H.

The analysis and discussion on the molecular signature of hydrophobicity of water at the interface with the oxides is conducted in Section 4, where the horizontal organization of interfacial water within a 2D-Hydrogen Bonded-Network is described and quantified from graphs.

The discussion on the oxide–water vertical organization of water is given in Section 5 through three sub-sections: two sub-sections (5.1 and 5.2) respectively reporting on aqueous SiO₂ and Al₂O₃, and sub-section 5.3 providing general considerations learned from these two aqueous oxides. The spectroscopic SFG signatures associated with these vertical motifs are discussed in Section 6. Section 7 provides a final discussion and opens perspectives.

2 Methodology

2.1 DFT-MD simulations

In this work, five silica oxide (SiO₂) and four alumina oxide (Al₂O₃) surfaces are modelled at the interface with liquid water, using DFT-based molecular dynamics simulations (DFT-MD or AIMD for Ab Initio MD). Among the SiO₂ surfaces, 3 of them are crystalline facets, 2 are amorphous facets with different surface silanol densities, while 2 Al₂O₃ surfaces are α-polymorphic and 2 are of the γ-polymorph. All surfaces are modelled as fully hydroxylated (i.e. terminated by hydroxyls, –OH), mimicking PZC (Point of Zero Charge) conditions. These conditions correspond to pH = 2–4 for silica oxides and pH = 6–9 for alumina oxides, where the net surface charge of the oxides is expected to be very close to neutral, making fully hydroxylated models of silica and alumina oxides reasonable representations.^37,44–48 Details of all DFT-MD simulation boxes are provided in Table 1 for both silica and alumina liquid–water interfaces. The oxide facets are reported by increasing the order of the silanol/aluminol surface densities. A scheme of the simulation cells is provided in Fig. 1, consisting of three elements: an oxide slab with the facet of interest oriented along the (xy) plane and ∼6–8 Å thickness in the vertical z-direction (∼2–3 layers of the solid), bulk liquid water of ∼20 Å thickness, and a vacuum layer above the water. This vacuum layer prevents the water from being confined between the periodic images of the solid part (i.e. the oxide) in the vertical direction (z-axis), and helps maintain the density of liquid water within the water slab. Two aqueous interfaces are hence simulated, one between the oxide top-surface and liquid water and one between the liquid water and the air, while the second oxide surface is in contact with the air in these simulation cells. Note that the simulation cells of the SiO₂ amorphous facets do not include the vacuum, two oxide–liquid water interfaces are hence simulated within the simulation boxes.

Fig. 1 Scheme of the simulation boxes employed in this work.

Table 1 Simulation cells of the DFT-MD trajectories conducted in this work for the SiO₂ and Al₂O₃ oxide facets at the interface with liquid water, presented by increasing value of the silanol/aluminol density at the top surface. Simulations done at PZC conditions (Point of Zero Charge). Top: SiO₂, bottom: Al₂O₃

Facet	Lateral surface area of the slab (nm^2)	Dimensions			Density of surface OH sites (number/nm^2)
		a (Å)	b (Å)	c (Å)
SiO2
am3.6	1.681	12.670	13.270	37.000	3.6
am4.7	1.490	9.117	16.342	32.000	4.7
(102)	2.081	14.897	13.967	45.000	5.8
(110)	2.276	15.766	14.439	45.000	7.9
(111)	2.439	14.100	17.300	45.100	9.8
Al2O3
α-(0001)	2.342	14.214	16.478	35.000	15.4
γ-(100)	2.787	16.723	16.663	55.000	15.4
γ-(110)	2.687	16.069	16.723	51.300	18.2
α-(10 1̄1)	3.459	19.036	18.172	60.000	18.3

---

The silica oxide–water interfaces include the (102), (110) and (111) α-crystalline SiO₂ facets, respectively of 5.8, 7.9 and 9.8 SiOH/nm² silanol densities, and amorphous facets of 3.6 and 4.7 SiOH/nm². The alumina facets are the α-(0001), the γ-(100), the γ-(110) and the α-(101̄1), with aluminol densities of respectively 15.4, 15.4, 18.2 and 18.3 AlOH/nm².

Fig. 2 and 3 present top and side views of the surfaces of selected silica (Fig. 2) and alumina (Fig. 3) oxides examined in this study. These representations highlight the diversity in both the spatial distribution and the chemical nature of surface silanol and aluminol groups at the oxide interfaces.

Fig. 2 Top (up) and side (bottom) views of two SiO₂/liquid–water interfaces studied in this work presented from the oxides top-surfaces only. Water molecules are excluded from the pictures. Geminal –Si(OH)₂ are depicted in green and isolated –Si(OH) in iceblue.

Fig. 3 Top (up) and side (bottom) views of two Al₂O₃/liquid–water interfaces studied in this work presented from the oxides top-surfaces only. Water molecules are excluded from the pictures. µ₁-OH₂ sites are depicted in orange, µ₁-OH in green, µ₂-OH in iceblue and µ₃-OH in cyan.

All DFT-MD simulations of the aqueous interfaces are carried out with the CP2K^49,50 software using GGA exchange–correlation functionals (BLYP ⁵¹ for the silica oxide–water systems, PBE ^52,53 for alumina oxide–water systems), augmented by Grimme D2 dispersion corrections.^54,55 We employ the Gaussian Plane-wave (GPW) scheme with GTH pseudopotentials⁵⁶ and a dual-basis set combining a plane-wave cutoff of 400 Ry (SiO₂) and 600 Ry (Al₂O₃) with Gaussian DZVP-MOLOPT-SR basis functions. Born–Oppenheimer molecular dynamics are performed at T = 300 K, using a 0.4 fs time step. Each system is equilibrated for ∼20 ps in the NVT ensemble (using the velocity rescaling scheme to ensure thermalization), followed by NVE production trajectories of ∼80 ps. These choices provide a good balance between accuracy and efficiency for the calculation of structural and vibrational properties at oxide–liquid water interfaces.^{29,37,47,48,57,58}

For the purpose of analyses and definition of the aqueous interfaces, the instantaneous interface is calculated at every MD step through the Willard and Chandler (WC) scheme.⁵⁹ Following our previous works,^29,37,48,60 the WC interface is used further to partition the water slab into the Binding Interfacial Layer (BIL), the Diffuse Layer (DL), and the bulk-like region. Following previous work,^29,37,48 joint probability distributions of H-bond lengths and donor O–H orientations (defined as the cosine between the O → H vector and the solid surface normal) are evaluated, allowing to extract BIL, DL and bulk water regions together with density values of the water. Bulk water indeed has H-bond distances within 2.65–2.80 Å with isotropic orientations, whereas interfacial layers exhibit directional ordering in the BIL and in the DL. Hydrogen bonds are identified using the Galli criterion: O⋯O distance ≤3.2 Å and ∠O–H⋯O in the 140–220° range.⁶¹

2.2 Vertical and horizontal orders of water at the interface with a solid: x_H descriptor of hydrophilicity/hydrophobicity

In ref. 15 we introduced the x_H descriptor to measure the degree of hydrophilicity or hydrophobicity of an aqueous interface, at the molecular level of representation. To that end, the H-bonding properties of the water molecules in the BIL (Binding Interfacial Layer) are used in order to partition the H-bonds into Horizontal (H) water–water H-bonds (i.e., oriented parallel to the surface) and Vertical (V) water–top-surface oxide H-bonds (regardless of water acting as donor or acceptor in the H-bond). See the scheme in Fig. 4 for illustration. At each time-step of the dynamics, the number of horizontal H-bonds (n_H) and of vertical H-bonds (n_V) are counted, and averaged over the whole trajectory.

Fig. 4 Scheme illustrating the interface between top surface solid oxide and liquid water in the BIL where Horizontal (red dashed lines) and Vertical (blue dashed lines) H-bonds can be formed.

The x_H molecular descriptor¹⁵ is then calculated as:

accounting for how much the water in the BIL is horizontally organized.

By definition, x_H ranges from 0 to 1, with 0 corresponding to a fully hydrophilic interface and 1 to a fully hydrophobic one. A value of x_H = 0 indeed requires a large n_V value, indicating an interface where numerous vertical (“hydrophilic”) surface–water H-bonds are formed. In contrast, values approaching 1 correspond to n_V = 0, meaning that water molecules preferentially H-bond with each other while avoiding H-binding sites at the surface, as expected for hydrophobic systems. A surface–water H-bond is included into n_V if it is long-lived for more than 70% of the simulation time-length, to ensure statistical relevance of the H-bond pattern. Note that the silanol/aluminol at the top-surface is taken as the anchor to the H-bond for the time-length criterion, allowing exchanges of water molecules in the BIL being H-bonded to it within the 70% time-length. In contrast, water–water H-bonds in the liquid and in the BIL are very short lived^36,58,62,63 and constantly exchange in between water molecules within the H-bonded network, thus the time-scale criterion is not applied in the n_H calculation.

2.3 Topological graphs at the atomistic level of representation and subgraphs of motifs

To characterize the horizontal and vertical H-bond networks formed between water molecules and water-top-surface oxide sites, respectively, graph-theoretical approaches are employed, following our previous works.^41,43,64 Within graph theory, a molecular system is represented as a mixed, colored graph (hereafter referred to as a 2D-MolGraph), where vertices correspond to atoms and edges correspond to interatomic interactions (covalent bonds, hydrogen bonds, ionic or organometallic interactions, etc.). In this work, we focus on the hydrogen-bond networks that can be formed within the Binding Interfacial Layer (BIL).^36,37,58 All interactions are represented by undirected edges, except for hydrogen bonds, which are represented by directed edges (arcs). Hydrogen atoms are not explicitly included as vertices of the 2D-MolGraphs. Instead, the direction of an H-bond edge encodes the identity of the heavy atom that is the donor of the H-bond: edges are directed from the donor heavy atom toward the acceptor heavy atom, as illustrated in Fig. 5. The hydrogen atom is therefore implicitly associated with the donor vertex.

Fig. 5 Schematic representation of a H-bond as a directed edge in a 2D-MolGraph illustrated for a water–water dimer. In the graph, the dashed red arrow is the directed edge (arc) associated to the H-bond and being oriented from donor to acceptor atoms. Hydrogen atoms are not included in the vertices of 2D-MolGraphs, the vertices associated with oxygen atoms are colored in red.

For each snapshot of a Molecular Dynamics (MD) trajectory, a 2D-MolGraph is constructed based on geometric criteria (interatomic distances and angles). For instance, a water–water H-bond is defined when the O⋯O distance is ≤3.2 Å and the ∠O–H⋯O lies between 140 and 220°.⁶¹ Changes in the molecular conformations sampled along the trajectory are quantified by monitoring changes in the graph topologies, which are tracked using graph-isomorphism algorithms.^41,65 Once all non-isomorphic 2D-MolGraphs have been identified over the whole trajectory, further graph analyses are performed on the set of non-isomorphic graphs in order to characterize the degree of collectivity of the water–water H-bonded network formed within the BIL on the one hand (i.e., horizontal organization of water in the BIL as described in Section 2.2), and in order to characterize the vertical organization of water in the BIL in terms of top-surface–water H-bonded motifs on the other hand (vertical organization of water in the BIL as described in Section 2.2).

Fig. 6a shows an example of a 2D-MolGraph of the BIL-water taken from one of our AIMD trajectories of an aqueous alumina oxide interface (for a given snapshot). In that graph, one can easily see that the water molecules (i.e. the red vertices in the 2D-MolGraph) are collectively H-bonded together, with only one water (one vertex) being isolated from the rest of the connected graph. All vertices, but one, are indeed interconnected to each other through directed edges that provide the direction of the H-bonds.

Fig. 6 (a) Illustration of a topological 2D-MolGraph for one snapshot extracted from a α-Al₂O₃/liquid water trajectory. Oxygen atoms are the vertices in red, hydrogen bonds are represented with the dashed red arcs (directed from the donor to the acceptor). (b) Distribution of the sizes of the connected components (CCs) of the ensemble of 2D-MolGraphs statistically analyzed from a MD trajectory (see text for details).

The degree of collectivity of the H-bonded network seen in this 2D-MolGraph can be quantified by the distribution of the size of the connected components (CCs) of the graph, i.e. the size of the subgraphs in which all vertices are mutually interconnected without any interruption. This is presented in Fig. 6b as the statistical analysis over the whole AIMD trajectory for the chosen Al₂O₃ aqueous interface presented in this illustration. The plot reports the percentage of appearance of each possible value for the size of the connected components within the graphs, the maximum size being the total number of water molecules located within the BIL. The distribution of the CCs for this particular example has one main peak centered at a CC size that corresponds to −85% of the number of water molecules that belong to the BIL. This is indeed the signature of a high degree of inter-connectivity between the water molecules within the BIL in contact with the oxide top-surface.

Beyond the measure of collectivity of the hydrogen bond networks formed between the water molecules in the BIL, one can characterize the vertical H-bonded network that can be formed between top-surface oxide silanols/aluminols and the water molecules in the BIL, directly from the 2D-MolGraphs. This will give us the knowledge of the structural patterns/motifs that exist in between top-surface silanols/aluminols and water, and their spatial distribution over the top-surface oxide. This motif analysis provides a geometrical mapping of the BIL and enables deeper insight into the structural, dynamical, spectroscopic, and reactive properties of aqueous interfaces. The knowledge of vertical motifs will ultimately be used in order to quantify the degree of hydrophilicity or hydrophobicity of the water wetting.

A motif (or 2D-Motif) is defined as a connected, induced subgraph of a 2D-MolGraph, characterized by its set of vertex types (chemical elements) and its internal bonding patterns (edges). Fig. 7c shows five representative motifs that could be found at the oxides/liquid water interfaces (the illustration is for Al₂O₃). At the top, one finds two simple motifs composed of one site of the surface acting as the donor of the H-bond (M311) or acceptor of the H-bond (M312) with one water molecule. Such motifs are denoted hereafter as 2-partner motifs composed of one top-surface aluminol (resp. silanol) and one water molecule. The sub-graphs are of size 3 (3 vertices). Below M311 and M312, one can then find in the figure two motifs of size 5, called 3-partner motifs, composed of 2 top-surface aluminol sites and one water molecule. The two surface sites are H-bonded to each other (follow the directed edge), and the water molecule can act either as a donor (M511) or acceptor (M512) in the H-bond with the aluminol (resp. silanol) top-surface site. The next motif reported below M511-M512, is a more elaborate motif between two H-bonded top-surface aluminols (resp. silanols) and two water molecules, one being donor and the other being acceptor of H-bonds with these surface sites. This is a 4-partner motif (M611) of size 6 (6 vertices).

Fig. 7 Illustration of the workflow from a 3D structure of top-surface oxide and BIL-water (left) to the associated 2D-MolGraph (middle) and to the extracted 2D-Motifs (right). The example is given for an α-alumina oxide/liquid water interface. In the 2D-MolGraph, vertices in brown are for aluminum atoms, vertices in red are for oxygens of water, vertices in orange are for oxygens of the top-surface alumina oxide. The sub-graphs of the 2D-MolGraph that are associated to pre-defined 2D-Motifs are shown on the right-side. Convention for the vertices are: brown for Al, orange for top-surface oxygen atoms, red for water oxygens. The direction of the edges in the 2D-Motifs provide the direction of the H-bond between water and top-surface.

These 2-, 3-, 4-partner motifs (2D-Motifs sub-graphs of respectively size 3, 5 and 6 with specific arcs) are the ones that will be recognized and quantified in this work for SiO₂– and Al₂O₃–liquid water BIL interfaces. Larger sizes of sub-graphs (size ≥7) as well as other combinations of directed edges (H-bonds) between the water molecules in larger motifs have not been investigated here.

A given 2D-Motif graph can appear multiple times within a given 2D-MolGraph, as the sub-graph can be associated to distinct combinations of vertices from the top-surface (but also from the water molecules in the BIL). Each appearance of a 2D-Motif is called an occurrence, regardless of the actual spatial position of this motif over the surface (i.e. regardless of the actual aluminol/silanol site at the top-surface that is partnered in the H-bond). Each 2D-Motif is characterized by its unique certificate, which is calculated based on the set of vertices. These certificates are reported in Fig. 7c with values M311-M611 for alumina oxides, and M301-M601 for silica oxides. The occurrences of a motif are found using the isomorphism algorithms.^41,65Note that the analysis of the 2D-Motifs has to be done with a hierarchical approach in order to remove any overlapping in the counting of the occurrences. Indeed, without this approach, a small motif (e.g. M311) that might be contained within a larger motif (e.g. M511) at the same location in space of the oxide aluminol/silanol sites, would be counted twice. However, the real motif to account for is the larger 511 one. Our analyses are therefore done by starting from the biggest size of motifs to identify (i.e., M611) and going down to the smallest size of motifs (e.g., M311). Whenever a small size motif is identified, a check is done in order to verify if this 2D-Motif has already been marked inside a larger size 2D-Motif. In the case of a positive check, the occurrence of the small size motif is not updated, otherwise the occurrence is updated and taken into account. All graph algorithms for the 2D-MolGraphs construction, isomorphism, and the 2D-Motif sub-graphs extraction are implemented in our GATEWAY ⁶⁶ and GRAPHETARIUM ⁴² home-made software.

2.4 SFG calculation

All details for SFG (Sum Frequency Generation) calculations from AIMD simulations have been given in our previous works, see for instance ref. 15,36,37,48 and 58. Here we just summarize the key elements.

The SFG signal is arising from the total resonant electric dipole non-linear susceptibility χ⁽²⁾(ω). Its real and imaginary components are calculated following the time-dependent method of Morita and Hynes,^67,68 using the model proposed by Khatib et al.⁶⁹ for dipole and polarisability derivatives. This model gives accurate SFG spectra both for band positions and absolute intensities, see e.g. some of our works.^{29,36,37,69,70} In more detail, based on the assumption that in the high frequency region (>3000 cm⁻¹) only the OH stretching motions are the main contributors to the signal, and neglecting intermolecular cross-correlation terms, the SFG signal is obtained by:

where (P, Q, R) are any x, y, z direction in the laboratory frame, and k_b and T are respectively the Boltzmann constant and temperature of the simulated system. 〈⋯〉 is a time-correlation function, α_PQ(t) and µ_R(0) are respectively the individual O–H bond contribution to the total polarization and dipole moment of the system and [small alpha, Greek, dot above]_PQ(t) and [small mu, Greek, dot above]_R(0) their time derivatives. M is the number of water molecules and n₁ and n₂ are two indices that identify each O–H oscillator per water molecule. In our works, ssp SFG signals are calculated, i.e. for the xxz directions. Using the direction cosine matrix (D) projecting the molecular frame (x, y, z) onto the laboratory frame (P, Q, R), and assuming that the O–H stretching is much faster than the modes involving a bond reorientation, one can write the two time derivatives as:

The D matrix and the projection of the velocities on the O–H bond axis (v_z) can be readily obtained from the DFT-MD trajectory, while (Raman tensor) and (APT) are parametrized from our previous works.^36,69,71,72

Hereby, the surfaces are neutral at PZC conditions, no SFG signal is therefore expected from DL layers (which we check in our calculations), the only contribution to SFG is arising from water in the BIL. Contributions to the SFG signal arising from top-surface hydroxyls are not included in this work.

3 What the x_H descriptor of hydrophilicity tells us about the aqueous interfaces of SiO₂ and Al₂O₃

Table 2 reports the x_H values of the molecular hydrophilicity descriptor calculated for the aqueous oxide interfaces investigated in this work. See Section 2.2 for the definition of x_H and its limiting values of 0 for hydrophilic and +1 for hydrophobic characters of the water in the Binding Interfacial Layer.

Table 2 Values of the x_H molecular descriptor of hydrophilicity for SiO₂ (top) and Al₂O₃ (bottom) liquid–water interfaces along with the density of hydroxyl groups (ρ_OH, in number of OH per nm²) at the top-surfaces. The presence or absence of the water horizontal hydrophobic 2D-HBonded-Network (2DN) is mentioned by symbols ‘✓’ or ‘–’, respectively. See Section 4 for the geometrical description of the horizontal hydrophobic 2DN

	SiO2
	am3.6	am4.7	(102)	(110)	(111)
xH	0.57	0.59	0.38	0.49	0.64
ρOH	3.6	4.7	5.8	7.9	9.8
2DN	✓	✓	—	✓	✓

---

	Al2O3
	α-(0001)	γ-(100)	γ-(110)	α-(10 1̄ 1)
xH	0.89	0.83	0.70	0.42
ρOH	15.4	15.4	18.2	18.3
2DN	✓	✓	—	—

---

Aqueous SiO₂ interfaces investigated here have x_H values varying between 0.38 and 0.64, i.e. from hydrophilic character to mildly hydrophobic. Only the (102) crystalline facet appears to have a clear hydrophilic wettability from the point of view of the x_H descriptor, all other facets (crystalline and amorphous alike) have a wettability character at the border between hydrophilic and hydrophobic (x_H = 0.49–0.64), and verging towards the hydrophobic side, regardless of the surface density ρ_OH in silanols.

On the other hand, three of the α- and γ-Al₂O₃ facets have a clear hydrophobic molecular character of wettability, with x_H = 0.89–0.70. Only α–(10 1̄ 1) is hydrophilic according to x_H.

These molecular hydrophobic characters of water at the interface with SiO₂ and Al₂O₃ oxides are somehow surprising when one considers that these oxides are characterized as macroscopically hydrophilic. It is also surprising to find a molecular hydrophobic wettability of these oxides if one considers the quite high surface density in silanols/aluminols at the top-surface of these oxides. If the SiO₂ amorphous facets have a low density of silanols at the surface (3.6 and 4.7 SiOH/nm²) that could explain that there cannot be too many vertical H-bonds being formed between the top-surface silanols and the water in the BIL, this is not obvious for the crystalline SiO₂ aqueous facets (ρ_OH = 6–10 SiOH/nm²), neither for the four Al₂O₃ aqueous facets which have large aluminol densities of 15–18 AlOH/nm². One would thus certainly have expected more hydrophilic behaviors (i.e., more vertical H-bonds) of these surfaces towards water.

There are of course other factors like the actual spatial distribution of silanols over the surface, their accessibility by water regarding the surface morphology and rugosity, and silanol/aluminol pK_a values that will affect the tendency of these sites to accept/donate H-bonds with water in the BIL. In previous works,^{10,47,57,73,74} pK_a values of silica and alumina top-surface hydroxyls were calculated, from AIMD simulations, showing in particular that aluminols are very basic species (pK_a = 16) in contrast to acidic silanols (pK_a = 5–10 with bimodal characters).

4 Horizontal order of water at the interface with oxides: a collective 2D-HBonded-network signature of molecular hydrophobicity

In previous works, we have shown^28,58,62 that hydrophobic water in the BIL is organized into a highly cooperative two-dimensional (2D) hydrogen-bonded network (2D-HBonded-Network, 2DN for short), with water–water hydrogen bonds being predominantly oriented parallel to the surface plane (and following its morphology). The emergence of this cooperative 2D-HBonded-Network constitutes a robust molecular signature of hydrophobicity that nicely correlates with macroscopic measurements of hydrophobicity.^15,28

The presence of the 2DN at any aqueous interface can be statistically analyzed by topological graphs, via the distribution of the connected components (CCs) within the graphs. See more discussions on the 2DN and CCs in topological graphs in our reviews^43,64,75 and in Section 2.3. The distributions were obtained from the statistical graph analysis of the DFT–MD trajectories. Plots of the distribution of the CCs of the topological graphs of the water in the BIL of silica and alumina oxides are reported in Fig. 8 and 9.

Fig. 8 Connected Components (CCs) from the statistical analysis of the 2D-MolGraphs of the water in the BIL of SiO₂ aqueous interfaces: (a) amorphous am3.6, (b) amorphous am4.7, (c) crystalline α-(102), (d) crystalline α-(110), (e) crystalline α-(111). The horizontal axis of these plots is the size of the CCs (the maximum value being the total number of water molecules in the BIL), the vertical axis is the percentage of each one of these CCs over the time-length of the AIMD trajectories.

Fig. 9 Connected Components (CCs) from the statistical analysis of the 2D-MolGraphs of the water in the BIL of Al₂O₃ aqueous interfaces: (a) α-(0001), (b) γ-(100), (c) γ-(110), (d) α-(10 1̄ 1). The horizontal axis of these plots is the size of the CCs (the maximum value being the total number of water molecules in the BIL), the vertical axis is the percentage of each one of these CCs over the time-length of the AIMD trajectories.

Whenever the water molecules in the BIL are organized as a 2DN, the distribution of CCs is composed of one main single peak located at a very large value of the connected components, that typically corresponds to more than 90% of the number of water molecules that belong to the BIL. This is the case of e.g. Fig. 8d for the SiO₂-(110) aqueous facet and of Fig. 9b for the γ-Al₂O₃-(100) aqueous facet. On the contrary, whenever the BIL-water is organized as hydrophilic, a rather homogeneous distribution of CCs is obtained from low to higher values of CCs, as can be seen in Fig. 8c (silica oxide) and 9c (alumina oxide).

Very interestingly, all the silica and alumina aqueous interfaces investigated in this work, except three of them – the SiO₂-(102) and the α-Al₂O₃-(10 1̄ 1) and γ-Al₂O₃-(110) facets – have their BIL-water organized into the collective 2D-H-Bonded-Network. Two of these aqueous interfaces have a x_H value of 0.38 and 0.42, that is clearly located on the hydrophilic wettability side of the x_H descriptor, while the BIL at the γ-Al₂O₃-(110) aqueous interface has a x_H value of 0.70, at the border between hydrophilic and hydrophobic wettability. The x_H descriptor nicely captures the balance between horizontal and vertical H-bonds made by the water molecules within the Binding Interfacial Layer. x_H values at the border between (clear) hydrophilic and (clear) hydrophobic characters, like the x_H = 0.70 for the γ-Al₂O₃-(110) aqueous interface, signifies that there are patches of hydrophobicity in the BIL but not to the extent of a very large inter-connected 2DN.

The existence of the water collective 2D-HBonded-Network in the BIL is hence the direct signature, in terms of intermolecular structural organization, of the molecular hydrophobic character of most of the silica and alumina oxide surfaces investigated in this work, when they are in contact with liquid water.

The collective 2D-HBonded-Network of the water in the BIL is composed of adjacent water–water H-bonded polygons of various sizes.^43,75 This can be observed in the topological graphs reported in Fig. 8 and 9 for one given snapshot taken out of the trajectories.

Fig. 10 directly illustrates this point for the BIL-water at the interface with the alumina oxide α-(0001) (Fig. 10a) and at the interface with the amorphous silica oxide am4.7 (Fig. 10b). In these figures, the statistical distribution of the size of the water–water H-bonded polygons within the 2DN is reported in percentages of occurrence. One can hence see that there is a predominance of water–water tetragons, pentagons and hexagons being formed within the collective 2DN network, that is the hallmark of the 2DN structural organization.^41,43,64,75

Fig. 10 Distribution in percentages of the size of water–water H-bonded polygons formed within the 2D-HBonded-Network: (a) Al₂O₃-α-(0001), (b) SiO₂ amorphous am4.7.

The distributions are centered at the size of the pentagons, with this latter having a higher probability of appearance. Note that we could not find a direct relationship between the oxide top-surface structural template made by the top-surface Si–O sites and the relative occurrence of these water–water H-bonded tetragons, pentagons and hexagons in the BIL-water.

5 Vertical order of water: surface–water motifs of hydrophilicity

While the water in the BIL at the interface with silica and alumina oxides can be organized into a collective hydrophobic 2D-HBonded-Network, as seen above, such organization does not preclude the water to form vertical H-bonds with the hydroxyls of the top-surface of these oxides. These vertical H-bonds are now analyzed with the topological graphs through the identification of pre-defined sub-graphs (2D-Motifs) of different sizes (in this work: 2-partner (sub-graph of size 3), 3-partner (sub-graph of size 5), and 4-partner (sub-graph of size 6) motifs as described in Section 2.3, in which, at least one partner belongs to the top-surface and another one belongs to water). Our aim is to find the ensemble of vertical top-surface silanols/aluminols–water H-bonded motifs that cover the oxide aqueous interfaces in the BIL at PZC conditions, to identify their complexity in size, to quantify their rate of appearance in time and space (100% or partial amounts of time? 100% or partial spatial covering of the top-surface?), to identify possible systematic recurrence of certain vertical motifs at oxide aqueous interfaces, and to identify whether certain geometrical criteria of the top-surface oxide could selectively induce the formation of specific vertical H-bonded motifs at the aqueous interface. Of course the pK_a chemistry of the silanols/aluminols at the surface is also influential in the formation of vertical H-bonds.

For each aqueous interface investigated in this work, we quantify vertical H-bonding motifs using two complementary approaches. First, we assess each motif appearance by the percentage of snapshots/2D-MolGraphs where each motif type appears at least once, regardless of which specific (Si/Al)–OH surface site participates in the motif (i.e., regardless of the specific location over the surface) or whether the participation of a given surface site changes between snapshots. A value of 100% indicates the presence of a given motif throughout the AIMD trajectory, at least once per snapshot. Second, we determine surface coverage in terms of the percentage of surface hydroxyl groups that are involved in each type of vertical motif. For this analysis, each Si–OH/Al–OH surface site is tracked individually (this is defined with a fixed ‘anchor labeling’) while the identity of the H-bonded water within the motif is allowed to vary over time. For each surface OH anchor, only occurrences of a given motif persisting for >70% (SiO₂)/>50% (Al₂O₃, for reasons discussed later in the text) are taken into account in the final statistics, hence ensuring the capture of stable HBonded motifs over a relevant time-length of the trajectory. We also remind the reader that a hierarchical analysis has been performed in the sub-graph analysis (see Section 2.3).

5.1 SiO₂ aqueous interfaces: co-existence of horizontal organization of water and its vertical anchoring

Fig. 11 (left) quantifies the percentage of appearance at least once per snapshot/2D-MolGraph of each vertical motif at SiO₂ aqueous interfaces. Simple 2-partner motifs (motifs 301 and 302, violet and orange bars) appear in 100% of snapshots across all interfaces, as expected for the most basic surface–water H-bonding pattern. The key finding is that these vertical silanol–water H-bonded patterns are appearing at all silica oxide aqueous interfaces, regardless of their hydrophilic or hydrophobic wettability character. While it seems natural to find such vertical H-bonding at the SiO₂-(102) hydrophilic facet, it is striking to find the same H-bonding patterns at all the hydrophobic facets. The larger size of top-surface–water vertical motifs are formed as well at the SiO₂ aqueous interfaces. One can hence see 3- (motifs 501 and 502, respectively in green and blue colors in the bars) and 4-partner motifs (motif 601 in red color in the bars) occurring at both hydrophilic and hydrophobic interfaces, such as at the crystalline (110) and (111) facets. At these facets, it is equally probable to find the water donating the H-bond in the 3-partner motif as receiving the H-bond. Interestingly, these 3- and 4-partner motifs can be seen with very low appearances (from 14% to 27%) at the hydrophilic SiO₂-(102) aqueous interface, while they have a 100% appearance at hydrophobic SiO₂ interfaces. Because of the low density of silanols at the surface of the amorphous SiO₂-3.6 SiOH/nm² facet, the inter-silanol distances exceed the H-bonding range, thus precluding the formation of 3- and 4-partner motifs.

Fig. 11 Analysis of the 2D-Motifs at the silica aqueous interface. Left: percentages of occurrence of a given 2D-Motif observed at least once per 2D-MolGraph (i.e. per snapshot in the AIMD trajectory). Right: percentages of Si–OH sites being involved in a vertical motif. Note that only the motifs that appear more than 70% of the DFT-MD/AIMD simulation time-length are taken into account. The color coding goes with the definition of the vertical motifs, here denoted by their certificate number. See Fig. 7 for the schemes of the motifs and definitions of the vertical H-bonded patterns.

The second analysis reported in Fig. 11 (right) shows the statistical percentages of Si-OH sites that are involved in the vertical 2D-Motifs at the SiO₂ aqueous interfaces. This analysis provides a quantitative view of the H-bonding coverage of top-surface hydroxyls by water in the BIL, and how many of the surface hydroxyls actually serve as anchors for vertical H-bonding, for both hydrophobic and hydrophilic aqueous interfaces. The percentage values are reported at the top of the colored bars.

While Fig. 11 (left) showed that simple 2-partner motifs appear at least once in 100% of snapshots across all interfaces, the detailed Si–OH site analysis in Fig. 11 (right) reveals striking heterogeneity in the actual surface coverage. The percentage of silanols forming persistent (i.e., for more than 70% of the trajectory duration) vertical H-bonded motifs varies from one interface to another: from complete involvement at am3.6 to only 59% at the geminal (111) facet. This nicely illustrates that motif appearance frequency does not equate to spatial coverage, highlighting the importance of distinguishing between instantaneous and persistent structural features. The low-density SiO₂-am3.6 surface exhibits complete vertical silanol coverage, while the high-density geminal crystalline SiO₂-(111) facet shows only 60% coverage silanols. This latter property is certainly due to over-crowded occupation of the surface by geminal silanols which leads to steric hindrance for the water molecules to approach and make H-bonds with both OHs per geminal site. Notably, silanols at non planar surfaces can simultaneously participate as both H-bond donors and acceptors, hence engaging simultaneously into two 2-partner motifs (motifs 301 and 302). This dual participation yields coverage exceeding 100% of the silanols as is the case at the amorphous am3.6 facet with a total percentage of 167% or with a total of 117% at the hydrophilic crystalline (102) facet, even though only 75% of the surface (102) silanols are actually involved in vertical motifs.

Counterintuitively, at the hydrophilic SiO₂-(102) aqueous interface, only simple 2-partner motifs are formed (with motif 301 dominating), despite compatible silanol–silanol spacing (2.70 Å) for inter-silanol H-bonding (mandatory for larger sizes of motifs). Amorphous hydrophobic interfaces (am3.6, am4.7) also display predominantly 2-partner motifs. At SiO₂-am4.7, the slightly higher silanol density and the local geometric organization of silanols into a patch however, enable 3-partner motifs (motif 501) covering 29% of surface sites. The 3- and 4-partner motifs (motifs 501, 502, 601) appear exclusively at crystalline (110) and (111) facets, requiring two conditions: sufficient silanol density with appropriate spacing (O–O distances of 2.32–3.09 Å at these surfaces) enabling inter-silanol H-bonds, and the organized horizontal 2D-HBonded-Network of water that provides horizontal H-bonded polygonal templates (i.e., water–water H-bonded tetragons, pentagons, hexagons) that spatially position water molecules for simultaneous horizontal (i.e., within the 2DN) and vertical bonding to the oxide surface. The hydrophobic 2DN thus appears to serve as a structural template that facilitates the formation of the larger vertical motifs. Finally, the 2DN anchors differently across surfaces: predominantly as a H-bond acceptor at SiO₂-am3.6 and (111), but as a H-bond donor at am4.7 and (110) silica oxide surfaces. This variability reflects both surface morphology and silanol acidity (pK_a values).

5.2 Al₂O₃ aqueous interfaces: vertical anchoring limited by pK_a surface chemistry

We applied the same analysis of vertical H-bonded motifs at the interface between Al₂O₃ and liquid water in the BIL. Fig. 12 (left) reports the percentage of occurrence at least once of the different motifs, and Fig. 12 (right) quantifies the percentage of involvement of the Al–OH sites in each vertical motif.

Fig. 12 Analysis of the 2D-Motifs at the alumina aqueous interface. Left: percentages of occurrence of a given 2D-Motif observed at least once per 2D-MolGraph (i.e. per snapshot in the AIMD trajectory). Right: percentages of Al–OH sites being involved in a vertical motif. Note that only the motifs that appear more than 50% of the DFT-MD/AIMD simulation time-length are taken into account. The color coding goes with the definition of the vertical motifs, here denoted by their certificate number. See Fig. 7 for the schemes of the motifs and definitions of the vertical H-bonded patterns.

Like the SiO₂, simple 2-partner motifs (motifs 311 and 312, violet and orange bars in Fig. 12, left) appear with 100% frequency at all Al₂O₃ interfaces, regardless of the hydrophilic (Al₂O₃-α-(10 1̄ 1)) or hydrophobic (all other interfaces) character of wettability. There is also a variability in the occurrences of the larger and more elaborated vertical motifs at the SiO₂ aqueous interfaces that does not exist at the Al₂O₃ aqueous interfaces. Moreover, while 3-partner motifs (motifs 511 and 512) appear instantaneously at all alumina oxide interfaces, no 4-partner motifs are formed. This indicates that although neighboring aluminols are close enough for their inter-surface H-bonding, the BIL water structure – whether hydrophobic or hydrophilic – cannot accommodate the spatial requirements for water molecules involved in 4-site vertical bonding.

The key distinction between alumina and silica oxide interfaces emerges when (time-)statistically examining the aluminols coverage by the vertical motifs. The higher basicity of top-surface aluminols¹⁰ results in weaker surface–water H-bonds and shorter lifetimes, requiring a relaxed persistence criterion of 50% (versus 70% for SiO₂) to capture persistent H-bonding features. Furthermore, alumina surfaces exhibit complex multilayered structures with µ₁-OHs/µ₁-H₂Os (top layer), µ₂-OHs/µ₂-O⁻s (sublayer), and µ₃-OHs/µ₃-O⁻s (deeper sublayer) in varying proportions (see ref. 48,76–82 for more details). The top-layered µ₁–OHs and µ₁–OH₂ are somehow largely (but not exclusively) dominating the capacity of the top-surface aluminols to form direct H-bonds with the water in the BIL. All together, these elements lead to 40–50% of aluminols that systematically remain uninvolved in vertical H-bonding across all Al₂O₃ aqueous interfaces.

The contrast between appearance frequency (Fig. 12, left) and persistent (Fig. 12, right) coverage is particularly striking for Al₂O₃. Ultimately, only simple 2-partner motifs achieve long-term stability at aqueous Al₂O₃ interfaces. While 3-partner motifs appeared in 100% of snapshots (Fig. 12, left), at least once, only 3–8% of aluminols participate in 3-partner configurations when applying the 50% time criterion (Fig. 12, right). 3-partner motifs form transiently at alumina interfaces but lack the stability observed at silica surfaces.

Unlike SiO₂, where hydrophobic surfaces show higher vertical coverage than hydrophilic ones, Al₂O₃ exhibits a milder trend. The hydrophilic α-(10 1̄ 1) interface achieves 60% aluminol coverage by vertical motifs, slightly surpassing the hydrophobic α-(0001) (50%) and γ-(100) (40%) interfaces. This difference suggests that the weaker H-bonds characteristic of surface aluminols, diminish the mechanistic coupling between the horizontal 2D-HBonded-Network formation and vertical anchoring. As with SiO₂, no universal donor/acceptor rule governs 2DN anchoring of Al₂O₃: the γ-(100) 2DN anchors predominantly as a H-bond acceptor (33% of Al–OH sites involved in motif 311), while α-(0001) shows equal donor/acceptor character of the 2DN vertical anchoring (22% of Al–OH sites involved in both motif 311 and motif 312).

5.3 Lessons from two aqueous oxides

Across the SiO₂ and Al₂O₃ aqueous interfaces investigated in this work, it is paradoxical and striking to observe that the hydrophilic facets – which one might intuitively expect to be predominantly connected to water through vertical hydrogen-bonded motifs at the molecular scale – are in fact involved in fewer vertical motifs than the hydrophobic aqueous facets. For example, the SiO₂-(102) facet exhibits a smaller fraction of top-surface silanols engaged in vertical motifs (75% in total) compared to the hydrophobic facets, which range from 100% of silanols involved at the amorphous am3.6 surface to approximately 60% at the geminal crystalline (111) facet.

This property highlights that the horizontal organization of water in the two-dimensional hydrogen-bonded network of hydrophobic aqueous interfaces is closely coupled to the formation of vertical oxide–water hydrogen bonds.

At the molecular scale, hydrophobicity hence does not mean that interfacial water in the BIL lacks hydrogen bonds with surface oxide sites. Rather, hydrophobicity emerges from a highly collective, horizontally organized hydrogen-bonded network of water in the BIL that is anchored to the oxide surface via a dense array of vertical surface–water hydrogen bonds.

Unsurprisingly, the more elaborate vertical 3- and 4-partner motifs (e.g., motifs 501, 502, and 601 at SiO₂ aqueous interfaces) are predominantly observed at crystalline facets that exhibit a higher density of top-surface hydroxyl groups, such as the SiO₂-(110) and (111) silica surfaces, and that also present favorable inter-site spacing between nearest hydroxyl neighbors. In addition, the acid–base chemistry (pK_a) of surface hydroxyl sites, as well as the nature of the hydroxyl groups present at the surface, play a significant role in the formation of these larger motifs. Consequently, because Al₂O₃ surface hydroxyls are more basic than SiO₂ silanols, vertical motifs of larger size are only weakly stable over time. Conversely, excessive hydroxyl crowding at geminal facets of SiO₂ (such as the (111) surface) appears to hinder the formation of these larger vertical motifs.

We further find that the existence of the more elaborate 3- and 4-partner vertical motifs benefits from a well-organized water structure within the Binding Interfacial Layer (BIL). Such a water–water templating environment is provided by the hydrophobic 2D-HBonded-Network (2DN), in which water molecules are hydrogen-bonded to one another through tetragonal, pentagonal, and hexagonal arrangements within the collective interfacial network. These water–water H-bonded polygons constitute effective spatial templates, positioning water molecules in configurations that are favorable for donating and/or accepting vertical hydrogen bonds with top-surface hydroxyl groups.

As a final note, across both SiO₂ and Al₂O₃ aqueous interfaces, no universal pattern was found to govern whether the hydrophobic 2DN anchors to these oxide surfaces as a H-bond donor or acceptor. At SiO₂, the 2DN anchors predominantly as an acceptor at am3.6 and (111), but as a donor at am4.7 and (110). Similarly, at Al₂O₃, the γ-(100) 2DN anchors primarily as an acceptor while α-(0001) shows balanced donor/acceptor character. This variability reflects surface-specific factors including morphology and hydroxyl acidity (pK_a values).

6 Spectroscopic SFG signatures of vertical motifs at the SiO₂ aqueous interfaces

In this Section, we present the calculated spectroscopic signatures of the vertical hydrogen-bonded motifs identified in Section 5.1 at the SiO₂ aqueous interfaces, expressed in terms of the Sum Frequency Generation (SFG) signal. SFG is one of the few vibrational spectroscopic techniques capable of probing water directly located at the interface with solids, as in the case of water within the BIL. Since SFG is insensitive to the horizontal organization of water, the calculations reported here account exclusively for the contributions from water molecules participating in the vertical motifs. In this Section, we focus exclusively on the SFG signatures at the aqueous SiO₂ interfaces. The SFG contributions are analyzed for each type of vertical motif, with the aim of potentially identifying each motif solely from its own SFG signature.

The HD-SFG spectrum, i.e. the imaginary part of the 2nd order susceptibility in the SSP polarization state, of the BIL of the five SiO₂–liquid water interfaces presented in this work, has been calculated from the DFT-MD trajectories. Only the water contributions are included in the SFG signal, see Section 2.4 for theoretical details. Results are presented in Fig. 13.The schemes indicate the water molecule dipole moment orientations, giving rise to either positive or negative signals. The normal to the surface is directed towards the surface, as also shown in the figure.

Fig. 13 Theoretical sum frequency generation (SFG) vibrational spectra () at PZC conditions for all SiO₂–liquid water interfaces investigated here. The water contribution is only taken into account in . Convention for , vector normal to the surface for SFG calculations, is also reported.

Only the water dipole moments that are oriented along the normal to the oxide surface contribute significantly to the SFG signatures in the SSP polarization, i.e. only the vertical motifs described in Section 5 contribute to SFG. The SFG theoretical model described in Section 2.4 is valid only to take into account the O–H stretching motion contributions (i.e., elongation and shortening of the O–H covalent bonds) of the water molecules in the BIL to the spectroscopic signals within the 3000–4000 cm⁻¹ domain.

For both crystalline and amorphous facets, the BIL-SFG signal is composed of one positive band at higher frequency (∼3500 cm⁻¹ for crystalline and ∼3600 cm⁻¹ for amorphous) and one negative band at lower frequency (∼3100–3300 cm⁻¹). As shown in ref. 28,37,58, the lower frequency band has a universal assignment due to water molecules oriented with their dipole moments outwards the normal to the surface and H-bonded to other water molecules (either located in the BIL or in the bulk liquid water). Considering this universal assignment, it is no surprise that the same H-bond strength and the same intensity of the band is obtained whatever the crystalline and/or amorphous SiO₂ facet. Note that some of these water molecules can also be involved in vertical motifs where they act as hydrogen bond acceptors from the surface silanols. However, the water giving rise to that negative band in the BIL of the amorphous am4.7 lead to a larger negative band intensity. Considering that SFG intensities are primarily correlated to alignment of the water OH groups along the normal to the surface, one can deduce that the water molecules in the BIL of the am4.7 might have a higher degree of alignment than at the other SiO₂ aqueous interfaces investigated in this work.

In contrast, the higher frequency positive SFG band, associated with water molecules pointing their dipoles towards the surface, reveals a bit more diversity in its central frequency and broadening. While this band is centered from ∼3400 to ∼3600 cm⁻¹ for the BIL-water of the three crystalline facets, it is blue-shifted to ∼3500–3650 cm⁻¹ for the BIL-water of the two amorphous facets. Once again, only the intensity of this band at the amorphous am4.7 differs strongly from the other facets, again possibly showing a stronger alignment of the water OH groups in the BIL of this interface along the normal to the surface. The assignment of this SFG positive band arises from the water molecules involved in vertical motifs in which water donates an H-bond to (in-plane) silanols at the SiO₂ surfaces. The variation in frequency of this H-bonded signature from ∼3400 to ∼3650 cm⁻¹ reveals that the strength of such vertical H-bond is non uniform among the SiO₂ aqueous interfaces. The amorphous facets lead to the weaker H-bond strengths, close to free OH signatures.

The subsequent discussions will reveal the SFG signatures of the vertical SiOH–water motifs depending on their 2-, 3- or 4-partner compositions. Fig. 14 and 15 report the spectroscopic deconvolution of the SFG signal in terms of contributions arising from the 2-partner motifs (motifs 301 and 302, Fig. 14) and 4-partner motifs (motif 601, Fig. 15).

Fig. 14 Theoretical SFG spectra of the 2-partner vertical motifs 301 (top) and 302 (bottom) at PZC conditions for all SiO₂–liquid water interfaces investigated here. The water contribution is only taken into account in . Schemes of the graphs of the vertical 2D-Motifs 301 and 302 are reported in the top left. Convention for , vector normal to the surface for SFG calculations, is also reported.

As expected, motif 301 (in which water in the BIL is accepting an H-bond from a donating surface silanol) contributes to the negative SFG band at ∼3200 cm⁻¹ for all SiO₂ facets. There are however two exceptions (crystalline (102) and amorphous am4.7) where a positive band at ∼3600 cm⁻¹ is also observed in the signal. There are two reasons for this unexpected signature. One reason arises from the difficulty to properly define the normal to the surface of a non planar amorphous surface, as SiO₂-am4.7. We have defined the vertical z-direction as the normal to the surface, which is incorrect for such a corrugated surface. Some angles are hence ill-defined, and can wrongly contribute to the positive band. Another reason observed at the SiO₂-(102) facet, is the flexibility in orientation over time of the water molecule within a given motif, flipping from an upward/outward direction over periods of time. While the motif is mostly recorded as motif 301 for the majority of time, hence contributing to the SFG negative band, it has some contributions that correspond to motif 302 over shorter periods of time, thus contributing to the SFG positive band, however recorded within motif 301 in the final deconvolution.

Regarding motif 302 motif (Fig. 14-bottom) in which water acts as a H-bond donor towards surface silanols, one can see that they are contributing to the positive SFG band for all crystalline facets. There is however a small negative (close to zero) contribution to the SFG signal of motifs 302 at the two amorphous surfaces, which can be due to the somehow ill-defined normal to the surface of these corrugated surfaces as well as to the flipping in the orientation of the water dipole over short periods of time.

The positive SFG bands arising from motif 302 are centered at around ∼3400 cm⁻¹, with the exception of the band obtained for am4.7, where a shift of almost 100 cm⁻¹ to the higher frequency domain (∼3500 cm⁻¹) can be seen. This indicates much weaker water–silanol H-bonds at this particular interface. One can also remark that there is a variety in the intensity of the positive band arising from motifs 302 among the five aqueous interfaces, that was not observed in the positive band of the final SFG spectrum. This reflects the fluctuation in the absolute number of motifs 302 being present at each interface, as well as the alignment of water’s dipole moments along the normal to the surface.

Fig. 15 now illustrates the SFG contributions arising from the four-partner motifs (motif 601), for the (110) and (111) silica–water interfaces (only two facets where this more elaborate motif can be found). In this motif, one water donates an H-bond to a silanol and another water molecule receives an H-bond from a silanol. From such a motif, one thus expects one positive (as in motif 302) and one negative (as in motif 301) contribution to the SFG, which is indeed seen in Fig. 15. Again, the actual intensity reflects the number of such motifs being present at the interfaces, i.e. 66% of the silanols at the interface with the SiO₂-(110) facet and only 4% of the silanols at the interface with the SiO₂-(111) facet. One is thus not surprised to see that the intensity of the two bands is much smaller for the (111) facet than for the (110) facet.

Fig. 15 Theoretical SFG spectra of the 4-partner vertical motif 601 at PZC conditions for the two SiO₂–liquid water interfaces that contain this motif. The water contribution is only taken into account in . Scheme of the graph of the vertical 2D-Motif 601 is reported in the top left. Schemes of orientation of the water molecules responsible for the bands are also presented. Convention for , vector normal to the surface for SFG calculations, is also reported.

The deconvolutions of the SFG signal in terms of the separate contributions of the different vertical motifs formed at the interface, also show the diversity of H-bonding patterns between water in the BIL and the silanols at the surface. Indeed, the deconvolved bands are rather broad, spanning over ∼100–150 cm⁻¹. However, in the end, water in both 2-partner and 4-partner motifs contribute to the same final SFG band position. It is thus hardly possible to tell whether a positive band in the SFG spectrum arises more likely from a 2-partner motif or from a 4-partner motif.

7 Final discussion on how water wets an oxide surface and associated SFG spectroscopic signatures

This work has given views on the molecular description of how water wets an oxide surface, taking SiO₂ silicon oxides and Al₂O₃ alumina oxides as two representative examples. Our work has focused on PZC conditions (Point of Zero Charge) of various crystalline and amorphous facets, for which the oxides' surfaces are composed of silanol (SiO₂) and aluminol (Al₂O₃) neutral sites.

Despite SiO₂ and Al₂O₃ being macroscopically known as hydrophilic surfaces, we have shown that the water in the Binding Interfacial Layer (BIL) of these two oxides is structurally organized as a highly collective 2D-HBonded-Network (2DN) with water–water H-bonds formed parallel to the oxide surface and following the morphology of the oxide surface. Such a collective H-bonded network of water at the interface with the solid is the molecular structural signature of hydrophobic water. Hydrophobicity of interfacial water was confirmed by the values of the x_H descriptor that measures the balance between the horizontal (i.e., hydrophobic water–water 2D-HBonded-Network) and vertical (i.e., hydrophilic top surface oxide–water H-bonds) organization of water in the BIL.

There are few exceptions to this hydrophobic organization of the water in the BIL within the series of SiO₂ and Al₂O₃ oxides facets (crystalline, amorphous, different polymorphs) investigated in this work. Only a third of these oxide surfaces were indeed found to be wetted by water with a definitive hydrophilic character at the microscopic level.

Strikingly, that water in the BIL is being organized as hydrophobic does not preclude the water to form vertical hydrophilic H-bonded motifs with the top surface oxide silanols/aluminols. We found that 100–60% (SiO₂ aqueous interfaces) and ∼60–40% (Al₂O₃ aqueous interfaces) of the surface silanols/aluminols are covered by these vertical H-bonded motifs.

There is therefore a balance being achieved between the horizontal water–water H-bonds and vertical top surface oxide–water H-bonds being formed in the BIL of SiO₂ and Al₂O₃ aqueous interfaces.

At hydrophobic aqueous interfaces, the 2D-HBonded-Network is hence anchored to the silanol/aluminol surface sites through at least 60% (SiO₂) of the top-surface sites. At hydrophilic SiO₂, where the 2D-Network is absent, the coverage in silanols by vertical H-bonded motifs is barely larger, achieving 75% of the top-surface silanols. There is therefore a dense coverage of silanols/aluminols by vertical H-bonded motifs at the aqueous oxides, whatever the actual wettability property of the surface. Note that the density in silanols/aluminols at the surface has no formal impact on the final coverage of the surface by the vertical motifs.

At the molecular level, the hydrophobic wettability of oxides manifests through the formation of a collective, horizontally organized 2D-HBonded-Network of water within the BIL. This network is anchored to the oxide surface via a dense array of vertically oriented hydrogen-bonded motifs. Such anchoring can arise from water molecules acting either predominantly as hydrogen-bond donors to the oxide surface or as hydrogen-bond acceptors, with no specific governing rule other than the acid–base (pK_a) chemistry of the surface hydroxyl groups. Conversely, the hydrophilic wettability of oxides at the molecular level is characterized by vertical hydrogen-bonded motifs that closely resemble those responsible for anchoring the 2D-HBonded-Network at hydrophobic aqueous interfaces.

Strikingly, we observe that the hydrophobic wetting of oxides controlled by the vertical anchoring of the 2DN is generally enabled by vertical motifs of greater size than those associated with hydrophilic wetting. We indeed observe that the water–water hydrogen-bonded polygons (primarily tetragons, pentagons, and hexagons) forming the hydrophobic collective 2D hydrogen-bonded network provide an effective structural basis for anchoring these vertical motifs. In particular, vertical motifs involving three or four partners – where two surface hydroxyl groups are hydrogen-bonded to each other and one or two water molecules are hydrogen-bonded to these hydroxyls – can be stabilized at the interface owing to the presence of these larger polygons within the horizontal water–water 2DN. These hydrogen-bonded polygons simultaneously organize neighboring water molecules geometrically and orient their O–H bonds in ways that favor the formation and stabilization of vertical anchoring between the oxide surface and the BIL water.

By contrast, hydrophilic wetting appears to rely predominantly on simpler vertical two-partner motifs composed of a single surface hydroxyl group hydrogen-bonded to a single water molecule.

As a final remark, we conclude that there is no inconsistency between macroscopic and microscopic wetting behaviors. Indeed, the SiO₂ and Al₂O₃ aqueous interfaces examined in this work are macroscopically hydrophilic, meaning that water readily wets these surfaces: this behavior is also observed microscopically, as we have shown that surface silanol and aluminol groups are densely covered by water molecules being hydrogen-bonded to these sites through vertical motifs. At the same time, the water in the Binding Interfacial Layer (BIL) of microscopically hydrophobic oxide aqueous interfaces exhibits a collective organization in which horizontal water–water hydrogen bonds are formed, leading to the formation of a 2D-HBonded-Network. This collective structuring of interfacial water constitutes a molecular signature of the hydrophobic character at oxide–water interfaces, directly captured by the molecular descriptor x_H: this 2DN however does not preclude these hydrophobic oxide surfaces to be wetted (H-bonded) vertically. This reconciles macroscopic and microscopic wettings.

Unfortunately, to the best of our knowledge, the current developments of Sum Frequency Generation (SFG) spectroscopy do not allow SFG to be sensitive to the horizontal organization of water within the interfacial 2DN. Instead, SFG primarily probes the vertical anchoring of water to oxide surface sites, providing information on both the strength of vertical hydrogen bonds and the orientations of water molecules within the vertical motifs.

The molecular manifestations of hydrophobic and hydrophilic wetting at oxide–water interfaces identified in this work naturally raise the question of their broader implications in physical and chemical processes that occur at oxide aqueous interfaces. These combined horizontal and vertical arrangements of interfacial water are expected to either facilitate or hinder some interfacial processes. For example, it is reasonable to anticipate that the transport of electrolytes, pollutants, or reactive species could be strongly influenced by these specific hydrogen-bonded interfacial structures. Similarly, catalytic and electrochemical reactions taking place at oxide–water interfaces are likely to be affected by such interfacial water organization. Investigations of these effects are already underway.^38,39

Our recently developed theoretical model for SFG spectroscopy, based on a database of vertical motifs and their populations – termed the “pop-model”⁸³ – can incorporate the SFG signatures of 2-, 3-, and 4-partner vertical H-bonded motifs identified in our studies, as presented in Section 6. The observed variability in the positions of the SFG bands for these motifs, which reflects both the dynamics and flexibility of water molecules within the BIL and the heterogeneity in the pK_a values of surface hydroxyls and surface morphologies, can be captured within the precompiled database of SFG signatures. This approach represents a step toward our ultimate goal: predicting both horizontal and vertical arrangements of interfacial water directly from the deconvolution of experimental SFG spectra using the pop-model, without the necessity of explicit simulations of aqueous interfaces.

Author contributions

Álvaro Cimas: methodology, writing – original draft, supervision, data curation, formal analysis, software, validation. Sana Bougueroua: methodology, writing – original draft, supervision, data curation, formal analysis, software, validation. Yijing Yang: investigation, formal analysis, data curation. Nargiz Damirova: investigation, formal analysis, data curation. Terry Chauvin: investigation, formal analysis, data curation. Evgenii Iablonovskii: investigation, formal analysis, data curation. Marie-Pierre Gaigeot: conceptualization, methodology, validation, writing – original draft, writing – review & editing, project administration, supervision, funding acquisition, resources.

Conflicts of interest

There are no conflicts to declare.

Data availability

Graph analyses have been performed with our homemade software GaTewAY, the description and code can be found at the following link: https://diamond-diadem.github.io/en/codes/scientific-computing/gateway/.⁶⁶ The description of our homemade Graphetarium software can be found in ref. 42. S. Bougueroua, A. Cimas, and M.-P. Gaigeot can also be directly contacted by e-mail to obtain access to graph software and data.

Acknowledgements

HPC resources from GENCI-France are acknowledged through the grant 072484 (CINES/IDRIS/TGCC). This work was partially funded by the France 2030 programme “ANR-11-IDEX-0003” by the Integrative Institute of Materials of Université Paris-Saclay – 2IM@UPSaclay. Prof Dominique Barth is acknowledged for all discussions regarding topological graphs and sub-graphs. Chloé Helain and Coralie Zens are acknowledged for initiating the work on sub-graphs. Dr Simone Pezzotti and Dr Wanlin Chen are acknowledged for their decisive contributions in the Gaigeot group in the early developments of the concepts of horizontal and vertical arrangements at aqueous interfaces.

Notes and references

1. J. L. Bañuelos, E. Borguet, G. E. J. Brown, R. T. Cygan, J. J. DeYoreo, P. M. Dove, M.-P. Gaigeot, F. M. Geiger, J. M. Gibbs, V. H. Grassian, A. G. Ilgen, Y.-S. Jun, N. Kabengi, L. Katz, J. D. Kubicki, J. Lützenkirchen, C. V. Putnis, R. C. Remsing, K. M. Rosso, G. Rother, M. Sulpizi, M. Villalobos and H. Zhang, Chem. Rev., 2023, 123, 6413–6544.
2. O. Björneholm, M. H. Hansen, A. Hodgson, L.-M. Liu, D. T. Limmer, A. Michaelides, P. Pedevilla, J. Rossmeisl, H. Shen, G. Tocci, E. Tyrode, M.-M. Walz, J. Werner and H. Bluhm, Chem. Rev., 2016, 116, 7698–7726.
3. D. J. Hoffman, S. W. Devlin, D. Garratt, S. Jamnuch, J. A. Spies, B. R. Nebgen, D. Schacher, A. Do, F. Bernal, E. J. Riffe, K. Kunnus, C. Y. Hampton, J. Duris, D. Cesar, N. Sudar, G. L. Dakovski, W. S. Drisdell, K. V. Lawler, A. Marinelli, M. W. Zuerch, R. J. Saykally, C. P. Schwartz, T. A. Pascal and J. D. Koralek, Nat. Comm., 2026, 16, 10522–10534.
4. J. Elsner, K. N. Lausch and J. Behler, Chem. Phys. Rev., 2025, 6, 031310.
5. Y.-B. Zhuang, C. Liu, J.-X. Zhu, J.-Y. Hu, J.-B. Le, J.-Q. Li, X.-J. Wen, X.-T. Fan, M. Jia, X.-Y. Li, A. Chen, L. Li, Z.-L. Lin, W.-H. Xu and J. Cheng, Sci Data, 2025, 12, 997–1096.
6. J. A. Nauman, D. Suvlu and A. P. Willard, Ann. Rev. Phys. Chem., 2025, 76, 181–202.
7. S. Sun, H. Yao, J. Pan and Z. Xian, J. Chem. Phys., 2024, 161, 214708.
8. C. Zhang, M. F. C. Andrade, Z. K. Goldsmith, A. S. Raman, Y. Li, P. M. Piaggi, X. Wu, R. Car and A. Selloni, Nat. Comm., 2024, 15, 10270.
9. A. P. Fellows, A. D. Duque, V. Balos, L. Lehmann, R. R. Netz, M. Wolf and M. Thamer, Langmuir, 2024, 40, 18760–18772.
10. M.-P. Gaigeot, M. Sprik and M. Sulpizi, J. Phys. Condens. Matter, 2012, 24, 124106.
11. X. Li, F. S. Brigiano, S. Pezzotti, X. Liu, W. Chen, H. Chen, Y. Li, H. Li, X. Lin, W. Zheng, Y. Wang, Y. R. Shen, M.-P. Gaigeot and W.-T. Liu, Nat. Chem., 2025, 17, 198–203.
12. F. Wei, S. Nihonyanagi and T. Tahara, J. Am. Chem. Soc., 2025, 147, 22402–22411.
13. L. Xinyi and L. Wei-Tao, Nat. Chem., 2025, 17, 167–168.
14. G. Gonella, E. H. G. Backus, Y. Nagata, D. J. Bonthuis, P. Loche, A. Schlaich, R. R. Netz, A. Kuhnle, I. T. McCrum, M. T. M. Koper, M. Wolf, B. Winter, G. Meijer, R. K. Campen and M. Bonn, Nat. Rev. Chem., 2021, 5, 466–485.
15. W. Chen, S. E. Sanders, B. Özdamar, D. Louaas, F. S. Brigiano, S. Pezzotti, P. B. Petersen and M.-P. Gaigeot, J. Phys. Chem. Lett., 2023, 14, 1301–1309.
16. Y. R. Shen, Second Harmonic and Sum-Frequency Spectroscopy, World Scientific, 2023.
17. Y. R. Shen, Fundamentals of Sum-Frequency Spectroscopy, Cambridge University Press, 2016.
18. C. Tian and Y. Shen, Surf. Sci. Rep., 2014, 69, 105–131.
19. Q. Du, R. Superfine, E. Freysz and Y. Shen, Phys. Rev. Lett., 1993, 70, 2313.
20. A. Eftekhari-Bafrooei and E. Borguet, J. Am. Chem. Soc., 2010, 132, 3756–3761.
21. Z. Zhang, L. Piatkowski, H. J. Bakker and M. Bonn, Nat. Chem., 2011, 3, 888–893.
22. C.-S. Hsieh, M. Okuno, J. Hunger, E. H. Backus, Y. Nagata and M. Bonn, Angew. Chem., Int. Ed., 2014, 53, 8146–8149.
23. A. M. Jubb, W. Hua and H. C. Allen, Annu. Rev. Phys. Chem., 2012, 63, 107–130.
24. L. Dalstein, E. Potapova and E. Tyrode, Phys. Chem. Chem. Phys., 2017, 19, 10343–10349.
25. S. Pezzotti, D. R. Galimberti and M.-P. Gaigeot, Phys. Chem. Chem. Phys., 2019, 21, 22188–22202.
26. F. Tang, T. Ohto, S. Sun, J. R. Rouxel, S. Imoto, E. H. Backus, S. Mukamel, M. Bonn and Y. Nagata, Chem. Rev., 2020, 120, 3633–3667.
27. Q. Du, E. Freysz and Y. R. Shen, Science, 1994, 264, 826–828.
28. S. Pezzotti, A. Serva, F. Sebastiani, F. S. Brigiano, D. R. Galimberti, L. Potier, S. Alfarano, G. Schwaab, M. Havenith and M.-P. Gaigeot, J. Phys. Chem. Lett., 2021, 12, 3827–3836.
29. J. D. Cyran, M. A. Donovan, D. Vollmer, F. S. Brigiano, S. Pezzotti, D. R. Galimberti, M.-P. Gaigeot, M. Bonn and E. H. Backus, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 1520–1525.
30. R. Godawat, S. N. Jamadagni and S. Garde, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 15119–15124.
31. A. J. Patel, P. Varilly and D. Chandler, J. Phys. Chem. B, 2010, 114, 1632–1637.
32. H. Acharya, S. Vembanur, S. N. Jamadagni and S. Garde, Faraday Discuss., 2010, 146, 353–365.
33. S. Gómez, N. Rojas-Valencia, S. A. Gómez, C. Cappelli, G. Merino and A. Restrepo, Chem. Sci., 2021, 12, 9233–9245.
34. S. Shin and A. P. Willard, J.Comput. Chem., 2018, 14, 461–465.
35. Y.-C. Wen, S. Zha, X. Liu, S. Yang, P. Guo, G. Shi, H. Fang, Y. R. Shen and C. Tian, Phys. Rev. Lett., 2016, 116, 016101.
36. S. Pezzotti, D. R. Galimberti, Y. R. Shen and M.-P. Gaigeot, Phys. Chem. Chem. Phys., 2018, 20, 5190–5199.
37. S. Pezzotti, D. R. Galimberti and M.-P. Gaigeot, J. Phys. Chem. Chem. Phys., 2019, 21, 22188–22202.
38. B. J. Mohammed and P. Simone, Chem. Sci., 2025, 16, 12823–12832.
39. M. B. Jassar, W. tao Liu and S. Pezzotti, 2025, preprint, arXiv:2506.18467, DOI:  DOI:10.48550/arXiv.2506.18467.
40. S. Pezzotti, W. Chen, F. Novelli, X. Yu, C. Hoberg and M. Havenith, Nat. Rev. Chem., 2025, 9, 481–494.
41. S. Bougueroua, R. Spezia, S. Pezzotti, S. Vial, F. Quessette, D. Barth and M.-P. Gaigeot, J. Chem. Phys., 2018, 149, 184102.
42. S. Bougueroua, A. A. Kolganov, C. Helain, C. Zens, D. Barth, E. A. Pidko and M.-P. Gaigeot, Phys. Chem. Chem. Phys., 2025, 27, 1298–1309.
43. S. Bougueroua, Y. Aboulfath, A. Cimas, A. Hashemi, A. E. Pidko, D. Barth and M.-P. Gaigeot, C. R. Chim., 2024, 27, 1–23.
44. J. N. Malin, J. G. Holland and F. M. Geiger, J. Phys. Chem. C, 2009, 113, 17795–17802.
45. Y. Duval, J. Mielczarski, O. Pokrovsky, E. Mielczarski and J. Ehrhardt, J. Phys. Chem. B, 2002, 106, 2937–2945.
46. P. M. Dove, Am. J. Sci., 1994, 294, 665–712.
47. M. Pfeiffer-Laplaud and M.-P. Gaigeot, J. Phys. Chem. C, 2016, 120, 14034–14047.
48. L. Potier, S. Pezzotti, M.-P. Gaigeot and A. Cimas, J. Phys. Chem. C, 2025, 129, 5990–6002.
49. J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing and J. Hutter, Comput. Phys. Commun., 2005, 167, 103–128.
50. J. Hutter, M. Iannuzzi, F. Schiffmann and J. VandeVondele, Wiley Interdiscip. Rev. Comput. Mol. Sci., 2014, 4, 15–25.
51. A. Becke, Opt. Phys., 1988, 38, 3098–3100.
52. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
53. Y. Zhang and W. Yang, Phys. Rev. Lett., 1998, 80, 890.
54. S. Grimme, J. Comput. Chem., 2004, 25, 1463–1473.
55. S. Grimme, J. Comput. Chem., 2006, 27, 1787–1799.
56. S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B Condens. Matter, 1996, 54, 1703.
57. M. Sulpizi, M.-P. Gaigeot and M. Sprik, J. Chem. Theor. Comput., 2012, 8, 1037–1047.
58. S. Pezzotti, D. R. Galimberti and M.-P. Gaigeot, J. Phys. Chem. Lett., 2017, 8, 3133–3141.
59. A. P. Willard and D. Chandler, J. Phys. Chem. B, 2010, 114, 1954–1958.
60. X. Li, F. S. Brigiano, S. Pezzotti, X. Liu, W. Chen, H. Chen, Y. Li, H. Li, X. Lin and W. Zheng, et al., J. Nat. Chem., 2025, 17, 198–203.
61. J. A. White, E. Schwegler, G. Galli and F. Gygi, J. Chem. Phys., 2000, 113, 4668–4673.
62. A. Serva, S. Pezzotti, S. Bougueroua, D. R. Galimberti and M.-P. Gaigeot, J. Mol. Struct., 2018, 1165, 71–78.
63. S. Pezzotti, A. Serva and M.-P. Gaigeot, J. Chem. Phys., 2018, 148, 174701.
64. S. Bougueroua, Y. Aboulfath, D. Barth and M.-P. Gaigeot, J. Mol. Phys., 2023, e2162456.
65. B. D. McKay and A. Piperno, J. Symb. Comput., 2014, 60, 94–112.
66. S. Bougueroua, F. Quessette, D. Barth and M.-P. Gaigeot, ChemRxiv, 2022, preprint,  DOI:10.26434/chemrxiv-2022-1d5x8.
67. A. Morita and J. T. Hynes, J. Phys. Chem. B, 2002, 106, 673–685.
68. A. Morita and T. Ishiyama, Phys. Chem. Chem. Phys., 2008, 10, 5801–5816.
69. R. Khatib, E. H. G. Backus, M. Bonn, M. J. Pérez-Haro, M.-P. Gaigeot and M. Sulpizi, Sci. Rep., 2016, 6, 24287.
70. S. Pezzotti, D. R. Galimberti, Y. R. Shen and M.-P. Gaigeot, Minerals, 2018, 8, 305–321.
71. S. A. Corcelli and J. L. Skinner, J. Phys. Chem. A, 2005, 109, 6154–6165.
72. F. Perakis, L. De Marco, A. Shalit, F. Tang, Z. R. Kann, T. D. Kühne, R. Torre, M. Bonn and Y. Nagata, Chem. Rev., 2016, 116, 7590–7607.
73. A. Cimas, F. Tielens, M. Sulpizi, M.-P. Gaigeot and D. Costa, J. Phys. Condens. Matter, 2014, 26, 244106.
74. M. Pfeiffer-Laplaud, M.-P. Gaigeot and M. Sulpizi, J. Phys. Chem. Lett., 2016, 7, 3229–3234.
75. R. AbouHaidar, S. Bougueroua, D. Duflot, M.-P. Gaigeot, B. Wyslouzil and C. Toubin, Farad. Disc., 2025, 258, 396–418.
76. L. Potier, PhD thesis, 2019, Université Paris-Saclay, NNT 2019SACLE032.
77. S. E. Mason, C. R. Iceman, T. P. Trainor and A. M. Chaka, Phys. Rev. B, 2010, 81, 125423.
78. P. Euzen, P. Raybaud, X. Krokidis, H. Toulhoat, J.-L. Le Loarer, J.-P. Jolivet and C. Froidefond, in Alumina, Wiley-VCH Verlag GmbH, 2008, pp. 1591–1677.
79. J. G. Catalano, C. Park, Z. Zhang and P. Fenter, Langmuir, 2006, 22, 4668–4673.
80. M. Digne, P. Sautet, P. Raybaud, P. Euzen and H. Toulhoat, J. Catal., 2002, 211, 1–5.
81. T. P. Trainor, P. J. Eng, G. E. Brown Jr, I. K. Robinson and M. De Santis, Surf. Sci., 2002, 496, 238–250.
82. P. J. Eng, T. P. Trainor, G. E. Brown Jr, G. A. Waychunas, M. Newville, S. R. Sutton and M. L. Rivers, Science, 2000, 288, 1029–1033.
83. W. Chen, D. Louaas, F. S. Brigiano, S. Pezzotti and M.-P. Gaigeot, J. Chem. Phys., 2024, 161, 144115.

---