2D nanoconfinement distorts the solvation structure of hydroxide but not of hydronium

Abstract

Understanding ion-specific behavior in nanoconfined water is essential for controlling charge transport and selectivity in two-dimensional membranes. Motivated by recent experiments revealing anomalous dielectric and transport behavior of water confined between hBN sheets, we use machine-learning-accelerated first-principles molecular dynamics to investigate the interfacial propensities of hydronium and hydroxide ions under similar confinement. We find that hydronium remains interfacial across all confinement regimes, whereas hydroxide shifts toward the interior as the environment becomes more bulk-like. This contrasting behavior reflects the combined influence of hydrogen bonding, interfacial water layering, and the polarization of hBN, which collectively stabilize hydronium at the surface while making hydroxide slightly more favorable within the structured interior. These findings expose an asymmetry in ion–surface coupling and establish a microscopic origin for hydronium's enhanced interfacial affinity. The results provide mechanistic insight into ion partitioning in two-dimensional channels and highlight the collective structuring of confined water as a key determinant of interfacial ion thermodynamics.

---

1 Introduction

Two-dimensional nanochannels have emerged as a promising platform for developing energy-efficient and environmentally safe membranes for molecular separation and filtration.^1–4 Fast water transport has been observed in graphitic nanoconduits, such as graphene nanochannels, motivating applications ranging from desalination^5,6 and nanofiltration^7,8 to energy harvesting⁹ and lab-on-a-chip¹⁰ technologies. Continued progress in carbon-based nanofluidics and other 2D materials requires precise characterization and mechanistic understanding of water flow enhancement and interfacial slip as functions of nanoscale confinement, surface curvature, and atomic-level surface chemistry.

Hexagonal boron nitride (hBN) has recently gained attention as a complementary 2D nanofluidic material due to its exceptional mechanical strength, chemical stability, and tunable band structure when integrated with graphene.¹¹ Like graphene, hBN features an atomically smooth, honeycomb lattice and exhibits high in-plane thermal conductivity¹² and resistance to oxidation,¹³ offering advantages in device stability, and biocompatibility.¹⁴ Despite these structural similarities, recent experiments and simulations have revealed unexpectedly large differences in aqueous transport: water experiences one to two orders of magnitude higher interfacial friction on hBN than on graphene, reflected in reduced slip lengths and higher friction coefficients.^15–20

These observations underscore the critical importance of molecular-scale interfacial behavior in 2D confinement. A detailed understanding of how ions interact with—and subsequently modify—the interface is essential, as ion-specific effects govern membrane selectivity,^21,22 energy-conversion efficiency,²³ and hydraulic resistance. Such unusual transport phenomena arise from the distinct physics of nanoconfinement, where water adopts layered arrangements rather than behaving as a bulk-like continuum,^24–26 and where ion hydration shells become distorted or partially dehydrated as they navigate narrow channels.^27–31

Recent experiments further reveal striking anomalies in water confined between hBN sheets, including an ultralow out-of-plane dielectric constant,³² a ferroelectric-like in-plane dielectric response, and proton conductivities approaching those of superionic liquids.³³ Yet the molecular mechanisms responsible for these behaviors remain poorly understood. Furthermore, hBN has traditionally been viewed as chemically inert and hydrophobic, but recent studies have shown it to acquire substantial negative surface charge at neutral pH, even in the absence of defects, due to the preferential adsorption of hydroxide ions.^34–37 This intrinsic charging mechanism stands in contrast to graphene, which remains largely neutral under comparable conditions, and suggests that the hBN–water interface hosts ion-specific chemistry governed by surface electronic structure rather than defects alone. Determining how hydronium and hydroxide distribute and reorganize near hBN interfaces is therefore central to understanding physics of confined aqueous solutions in hBN-based systems.

In light of these advances, it is clear that ion specificity at aqueous interfaces arises from a delicate interplay among hydrogen bond topology, interfacial polarization, and confinement-induced solvent structuring. Yet, despite extensive work on graphene³⁸ and the air–water interface,^39–41 the microscopic behavior of hydronium and hydroxide at hBN–water interfaces under 2D confinement remains largely unexplored. This gap persists even as recent studies have begun to clarify other ion interactions and polarization effects at hBN interfaces.^42–44 In this work, we conduct a comprehensive theoretical study to elucidate the interfacial propensity of hydronium and hydroxide within nanoconfined water adjacent to hBN. Leveraging machine-learning-accelerated molecular dynamics trained on ab initio data, we achieve first-principles fidelity while accessing simulation scales and sampling depths far beyond traditional electronic structure methods. By quantifying ion distribution, solvation structure, hydrogen bond reorganization, and electron density response across confinement regimes, we provide molecular level insights into the origin of ion-specific stabilization and surface charging in hBN nanochannels. These results establish a mechanistic foundation for understanding proton and hydroxide behavior in 2D confined aqueous systems and inform the design of next-generation hBN based membranes and nanofluidic devices.

2 Methods

2.1 System set up

The systems studied here have two parallel free-standing hBN sheets, each containing 60 atoms, separated by heights ranging from approximately 6.6 Å to 20 Å. These different slit widths correspond to varying amounts of water molecules corresponding to system sizes between 17 and 94, all determined to give a water density of 1 g cm⁻³ between the nanoslits. The water-filling procedure used here follows Dufils et al.,⁴⁵ which reports dielectric responses in close agreement with experiment. All systems were simulated in orthorhombic simulation cells employing periodic boundary conditions in all three directions. For interfacial systems, we used a vacuum of 15 Å to uncouple periodic images in the z direction, leading to negligible interactions between the images. The lengths of the box in the x and y directions for all systems were 13.02 Å and 12.53 Å, respectively.

2.2 Machine learning potentials

In this work, we use the MACE architecture,⁴⁶ which allows for efficient data-efficient training with high-order equivariant message passing and has been proven robust in a wide variety of scenarios. We developed and validated a MACE machine learning potential (MLP) model with two message passing layers, a 6 Å cutoff distance, and 128 equivariant messages with body order 4 and spherical harmonics of maximal degree 1. The MLP captures semilocal interactions through a receptive field that spans the product of the number of layers and the cutoff distance per layer. In this case, the total receptive field is 12 Å, allowing the MLP to account for interactions within this range. While the model does not explicitly account for long-range effects, the 12 Å receptive field spans nearly the entire width of the slit in most cases, effectively capturing the relevant electrostatic interactions within the simulation. The final energy and force validation root-mean-square errors were 7 meV per atom and 100 meV Å⁻¹, respectively.

To accurately represent the potential energy surface of the systems, we train our MLP model using energies and atomic forces obtained from DFT calculations using the CP2K/Quickstep code.⁴⁷ We specifically used the revPBE0-D3 (ref. 48 and 49) functional due to its robust performance in reproducing the structure and dynamics of liquid water,⁵⁰ while also accurately capturing interactions between water and hBN surfaces.⁴⁵ Atomic cores are represented using dual-space GTH pseudopotentials.⁵¹ The Kohn–Sham orbitals of oxygen and hydrogen atoms are expanded using the TZV2P basis set, while those of boron and nitrogen atoms are expanded using the DZVP basis set. An auxiliary plane-wave basis with a cutoff of 750 Ry was used to represent the density. Training data was generated by first equilibrating systems for 400 ps with slit widths of 6, 10, 20, and 40 Å using the MACE foundational model,⁵² and then sampling from short 10 ps ab initio molecular dynamics (AIMD) runs with CP2K under the constant volume canonical ensemble (NVT) where the temperature is controlled by stochastic velocity rescaling⁵³ with a relaxation time of 1 ps and a time step of 0.5 fs. The final model consisted of 1931 structures, with 1572 involving hBN interfaces and 359 associated with bulk conditions. The configurations included are neutral, hydronium, and hydroxide systems both in bulk water and under confinement at slit widths of 6, 10, 20, and 40 Å For the protonic defect systems, the system's net charge was neutralized using a uniformly charged background. Parity plots of forces for configurations excluded from the training set, along with RDFs compared to short 10 ps AIMD simulations, are provided in the SI (Fig. S1–S3).

2.3 Molecular dynamics simulations

All MD simulations reported herein, which were based on the MLP, were performed using the ASE software⁵⁴ at a temperature of 300 K, in the NVT ensemble. Classical equations of motion were numerically integrated using the velocity Verlet algorithm⁵⁵ with a 0.5 fs time step. Temperature was controlled by a white-noise Langevin thermostat with a friction parameter of 0.01 fs⁻¹. Simulations were carried out for four slit widths: 6.63, 10.55, 15.55, 20.70 Å. For each of the four slit widths and two species, we conducted three independent simulations. Each simulation included a 100 ps equilibration period followed by a 3 ns production run. Uncertainties in reported values were calculated using the standard deviation from these replicates. All systems were simulated in orthorhombic simulation cells employing periodic boundary conditions in all three directions. The simulation cells were initially set up by randomly packing water molecules in each of the slit widths with a density of 1 g cm⁻³. To prevent interactions between the periodic images, 15 Å vacuum (exceeding the model's receptive field) was added in the z direction of these initial configurations. The hBN walls are kept fixed in the z direction (with thermal fluctuations allowed in the xy plane), following previous confined-water hBN studies^45,56 and because free-standing hBN shows only minimal intrinsic rippling (1–2 Å),¹⁷ making a fixed sheet a physically reasonable choice.

2.4 Electronic structure calculations

The electronic properties of the protonic defects at the hBN–water interface were analyzed using the same electronic structure settings used to train our MLP. To reduce the computational cost, revPBE⁴⁸ was used as the exchange correlation functional. To assess the interactions between the liquid environments and the hBN layers, we used DFT to analyze their electron density difference (Δρ), defined as Δρ = ρ_liq/hBN − ρ_hBN − ρ_liq, where ρ_liq/hBN, ρ_hBN, and ρ_liq are the electron densities of the system under consideration, the isolated hBN surfaces, and the isolated liquid environment, respectively. Atomic charge distributions were calculated using the Hirshfeld⁵⁷ charge analysis with CP2K.⁴⁷

3 Results and discussion

3.1 Interfacial propensity of water's ions under hBN confinement

Fig. 1 shows the oxygen density profiles for neutral confined water and for systems containing either a hydronium or hydroxide ion in aqueous solution. At the smallest slit width, the hydronium ion exhibits a clear preference for the interface, with its density significantly peaked near the hBN surface. While hydroxide, at the same narrow confinement, predominantly occupies the center of the water slab, avoiding the interface and instead localizing within the most bulk-like region available. As the confinement distance increases, this behavior somewhat persists. Hydronium consistently remains interfacial across all sampled slit widths, showing minimal population in the intermediate or interior water layers. Hydroxide, however, increasingly samples both the interface and the interior as confinement relaxes. While hydroxide maintains significant affinity for the surface, it also exhibits a pronounced probability of residing within the ordered, layered water structure that develops away from the surfaces. Interestingly, this trend parallels previous observations for water confined between graphene sheets, where hydroxide also demonstrated greater interior affinity than hydronium.³⁸ However, the effect appears reduced in hBN nanoconfinement. The hydroxide ion shows a somewhat weaker drive toward the interior compared to graphene systems, suggesting subtle differences in ion–surface interactions between the two 2D materials, potentially arising from hBN's partial polarity and the associated interfacial electrostatic environment. This observation is in line with predictions from Wang et al.³⁴ of a strong affinity of both hydroxide and hydronium ions for a water hBN surface, which they confirm experimentally to give rise to a negatively (positively) charged interface under mildly basic (acidic) conditions.

Fig. 1 Probability density profile of oxygen atoms in the system as a function of the z position for each sampled confinement distance. The red line represents the density profile for the neutral water system while the yellow represents the density profile of the hydronium ion and the green represents the density profile of the hydroxide ion. Z = 0 represents the midpoint between the two hBN surfaces. The dashed blue line represents the average positions of the nitrogen and boron atoms. Shaded regions represent the standard deviation over the three runs.

This interfacial asymmetry between hydronium and hydroxide echoes—though does not precisely replicate—their behavior at other aqueous interfaces. In graphene nanochannels, prior work has shown that hydronium strongly localizes at the surface, whereas hydroxide preferentially occupies interior layers, reflecting a confinement-amplified expression of their intrinsic solvation asymmetry.³⁸ The trends observed here for hBN are qualitatively similar but exhibit important distinctions: hydroxide shows a somewhat stronger interfacial presence than in graphene, consistent with reports that defect-free hBN can acquire negative charge through hydroxide adsorption. We expect that these differences arise from the distinct electronic and hydrogen-bonding environments of the two materials. Graphene presents an electronically homogeneous, nonpolar surface that disfavors hydroxide's anisotropic solvation motif, while hBN introduces atomic-scale polarity capable of stabilizing charge redistribution, thereby reducing the energetic cost of hydroxide approaching the interface. A related depth-dependent ion organization is also found at the air–water interface, where the acid–base character of the free surface has been the subject of longstanding debate.^39–41 Macroscopic charge-sensing experiments historically pointed to hydroxide enrichment, whereas molecular spectroscopies often indicated hydronium accumulation. More recent deep-potential molecular dynamics simulations⁴¹ add further nuance as they find a double-layer structure in which both ions accumulate near the interface, but at distinct depths. Taken together, these observations point to a unifying picture across interfaces of increasing complexity: hydronium consistently stabilizes at the liquid boundary, while hydroxide prefers regions where water maintains a more bulk-like coordination environment. The degree to which hydroxide can approach the interface is modulated by the local polarity and electronic structure of the confining surface. This consistency across air, graphene, and hBN interfaces underscores the fundamental role of hydrogen-bond topology and interfacial polarization in shaping ion propensities under confinement.

3.2 Asymmetric hydrogen-bond environments under nanoconfinement

To further elucidate the differing interfacial propensities of hydronium and hydroxide, we first examine how confinement alters their local hydrogen-bonding environments. Fig. 2 reports the average number of hydrogen bonds donated and accepted by each ion as a function of slit width, with bulk water values included for reference. The water regions corresponding to distinct layering are defined based on the density profiles shown in Fig. S4. Overall, the solvation structure of hydronium is remarkably robust to confinement. Across all sampled slit widths, hydronium maintains nearly the same number of hydrogen bonds as in bulk water, with a small reduction in accepted bonds at the tightest confinement. This indicates that even when strongly confined or positioned directly at the hBN surface, hydronium can largely preserve its characteristic three-donor hydrogen bond structure. In contrast, hydroxide exhibits pronounced confinement-induced restructuring. At the smallest slit width, hydroxide loses a substantial fraction of its accepted hydrogen bonds and essentially ceases to donate hydrogen bonds altogether. With increasing confinement distance, the hydrogen-bond coordination gradually approaches the bulk value. At the largest slit width, hydroxide recovers its bulk-like hydrogen-bonding behavior in the central “bulk-like” region of the water slab. However, when hydroxide resides at the interface at this widest spacing, it still shows a loss of hydrogen-bond donation, indicating a persistent asymmetry in solvation near the surface.

Fig. 2 Average number of hydrogen bonds donated and accepted by (top) hydronium and (bottom) hydroxide ions as a function of the four sampled confinement distances. The dashed lines represent the respective bulk water values. For the two larger confinement distances the number of hydrogen bonds for the intermediate and bulk-like regions are estimated as well. Hydrogen bonds are defined using the Luzar and Chandler definition.⁵⁸

These observations establish a fundamental asymmetry in how the two ions interact with the hBN interface under confinement. Hydronium can adsorb to the interface while retaining nearly bulk-like hydrogen-bonding. Hydroxide, by contrast, experiences significant disruption to its preferred hydrogen-bonding environment when positioned near the surface; this structural frustration leads to a decrease in interfacial localization, particularly at narrow confinement (e.g. 6.6 Å). Thus, the hydrogen-bonding analysis provides a microscopic explanation for the higher interfacial affinity of hydronium relative to hydroxide. Hydronium's ability to preserve its solvation network enables favorable adsorption without sacrificing hydrogen-bonding stability, whereas hydroxide's inability to maintain its solvation structure drives it towards the more solvated interior of the nanochannel. A similar asymmetry between water ions has been reported in water confined between graphene sheets,³⁸ reinforcing the generality of this behavior in two-dimensional nanoconfinement environments.

3.3 Orientation and site-specific interactions of water ions at hBN interfaces

To gain molecular-level insight into the surface-specific interactions governing ion stabilization at the hBN interface, we analyzed the radial distribution functions (RDFs) between the oxygen atom of each ion and the boron and nitrogen atoms of the surface (Fig. 3). Across all confinement distances, the hydronium ion consistently resides closer to the interface, with its oxygen atom located primarily within 2.5–3.0 Å of the surface, whereas the hydroxide ion remains slightly farther away, typically beyond 3 Å. This trend aligns with the density profiles in Fig. 1, which also show that hydronium maintains a stronger interfacial localization than hydroxide.

Fig. 3 Radial distribution functions of the oxygen atom of the respective ions with the boron and nitrogen atoms of the hBN surfaces (only accounted for when the ion is present at the interface of the system). Gray dashed lines are included as visual guides to highlight that the hydronium RDF peaks occur at shorter distances. Shaded regions represent the standard deviation over the three runs.

The RDFs further reveal distinct coordination characteristics for the two surface species. Although the oxygen–boron correlation peaks at slightly shorter separations, the first peak corresponding to the oxygen–nitrogen pair exhibits greater intensity for both ions, indicating a more favorable short-range interaction with nitrogen sites. In contrast, the oxygen–boron RDF displays a weaker first peak but a pronounced second maximum, suggesting secondary coordination at larger distances. When the ions are interfacial, these structural motifs remain largely invariant with confinement, and at the widest slit width the RDFs of hydronium and hydroxide become nearly indistinguishable. Taken together, these results demonstrate that both ions, particularly hydronium, exhibit preferential stabilization through close-range interactions with the nitrogen atoms of the hBN surface, reflecting the subtle polarity and electronic heterogeneity inherent to this two-dimensional material.

For both ions, the first O–N and O–B RDF peaks remain essentially unchanged with confinement, indicating that the preferred nearest-neighbor distances to the hBN surface atoms are stable whenever the ion is interfacial. However, the broader structural features evolve differently for the two water ions. Hydronium retains a strong interfacial preference across all slit widths, and at the smallest confinements its slightly increased separation from the surface leads to broader, less ordered second-shell structure in both the O–N and O–B RDFs. In contrast, although hydroxide's interfacial preference varies with confinement, its RDFs show more sharply defined second-shell features when it is at the interface, reflecting a more ordered positional relationship with the surface atoms under the same conditions. These differences are significant because they highlight that hydronium and hydroxide stabilize at the interface through distinct structural mechanisms, which contribute to their contrasting interfacial propensities under confinement.

Fig. 4 shows the orientation distributions of the O–H bond vectors with respect to the surface normal for the smallest and largest confinement distances. For hydronium, when located at the interface, all O–H bonds lie predominantly parallel to the hBN surface, producing a sharp peak near cos(θ) = 0 under 6.6 Å confinement. Then under 20 Å confinement there is a peak at cos(θ) = 0 and plateau from cos(θ) = −1 to 0. This indicates that under 20 Å confinement the hydronium O–H bonds at the interface either point parallel to the hBN surface or away from either surface. Hydroxide, in contrast, displays a broader and more complex angular distribution: although under 6.6 Å confinement it also shows a peak near cos(θ) = 0, it exhibits pronounced maxima at cos(θ) = ±1. While under 20 Å confinement the O–H bond has a pronounced peak at cos(θ) = 1 with a slight plateau from cos(θ) = −1 to 0. This pattern indicates that the hydroxide O–H bond frequently aligns along the surface normal, pointing towards either hBN wall. This orientation reflects the hydrogen-bonding asymmetry of hydroxide. Its donated O–H bond is relatively weak, making it easier to break than the stronger hydrogen bonds that hydroxide accepts, yet the ion simultaneously prefers a roughly square-planar solvation structure. Enforcing such a geometry in tightly confined, quasi-2D water is energetically costly. As confinement widens (20 Å), interfacial water becomes slightly more bulk-like, reducing this penalty and allowing the cos(θ) = ±1 orientations to appear more frequently. Overall, these trends suggest that hydroxide experiences a slightly larger energetic cost at the interface because its preferred orientation and solvation structure are more difficult to accommodate, whereas hydronium can reorient to maintain its hydrogen-bonding network with relatively little penalty.

Fig. 4 Cos(θ) angle of the respective OH bonds for either the hydroxide or hydronium ion with respect to the normal direction of the hBN surface. For 6.6 Å orientations are calculated with respect to the upper hBN surface. For 20 Å orientations are calculated with respect to either the upper or lower hBN surface depending on if the ion is located below or above the midpoint of the system. The definition of the angle θ is shown on the respective hydronium and hydroxide snapshots. Shaded regions represent the standard deviation over the three runs. For the 20 Å case the three different lines plotted represent when the ion is at the interface, or present in the intermediate region, or present in the bulk-like region.

3.4 Interfacial polarization drives ion adsorption at hBN

Despite the hydroxide ion exhibiting an undersolvated hydrogen-bonding environment at the interface, it appears that both hydronium and hydroxide ions are attracted to the hBN surface, as shown by Fig. 1 and 3. To elucidate the microscopic origin of this interfacial stabilization, electron-density difference profiles were computed for the 6 Å and 15 Å confinement distances. These calculations isolate the redistribution of electronic charge induced by the presence of the liquid and the solvated ion near the hBN surface, enabling direct visualization of how each species perturbs the electronic structure of the interface. Fig. 5 presents the integrated electron-density differences across the z direction at confinement distances of 6.6 and 15 Å, for hydronium and hydroxide, respectively. Notably, the qualitative behavior is consistent across both confinement regimes, indicating that the electronic response at the interface is robust with respect to separation.

Fig. 5 Constant electron density differences to quantify the interaction of hBN with water. (a) Isosurfaces representing constant electron density differences for a representative snapshot of the 15 Å hydroxide system. Blue isosurface indicates a region of electron depletion, whereas yellow isosurfaces indicate regions of electron accumulation. (b) Electron density difference profiles (integrated over the XY coordinates) for the 6 Å and 15 Å case for the hydronium and hydroxide system. The electron density difference is defined as Δρ = ρ_liq/hBN − ρ_hBN − ρ_liq, where ρ_liq/hBN, ρ_hBN, and ρ_liq are the electron densities of the system under consideration, the isolated hBN surfaces, and the isolated liquid environment, respectively. Profiles are averaged over 500 frames randomly sampled from the entire trajectory. Dashed lines represent constant electron density difference profiles while continuous lines represent oxygen atom density profiles (calculated from the 500 frames). The vertical blue line represents the position of hBN surface.

In both cases, the interaction between the aqueous environment and the hBN surface generates a clear polarization signature. For hydronium, the electron-density profile reveals a region of enhanced electron accumulation immediately adjacent to the surface, followed by a depletion layer farther into the liquid. In contrast, hydroxide exhibits the opposite pattern—an initial electron depletion at the interface followed by an accumulation region. This opposite-phase polarization response reflects the distinct charge distributions and solvation structures of the two ions, and highlights that the interface accommodates each ion through complementary electronic rearrangements. This analysis reveals pronounced image charge interactions, with the highly polarizable hBN surface generating induced electrostatic responses that attract ions of either sign. As shown in Fig. S6, the Hirshfeld charge distributions of nitrogen and boron atoms broaden substantially when the hBN sheet interacts with confined water with the ions, with boron exhibiting the most pronounced widening. Collectively, these results indicate that both ions elicit favorable polarization at the water–hBN interface, accounting for their persistent interfacial preference across all confinement distances. This interfacial stabilization, despite hydroxide's reduced hydrogen-bond coordination relative to bulk solvation, underscores the important role of surface-mediated electronic responses in governing ion adsorption in two-dimensional nanoconfinement.

3.5 Structural origins of ion spatial propensity under confinement

To better understand why hydroxide exhibits some propensity for the inner water layers at larger confinement distances, whereas hydronium shows no such preference, we examine how each ion perturbs the surrounding water structure. Fig. 6 shows the probability density distributions of the tetrahedral order parameter for water molecules⁵⁹ in the bulk-like region at a confinement distance of 20 Å. The distributions are compared across three cases, when the ion resides in the interfacial, intermediate, or bulk-like region, and are plotted alongside the corresponding distribution for the neutral water system. The corresponding probability density curves for water molecules in the interfacial and intermediate regions are provided in the SI (Fig. S5) and reveal no significant structural differences in the surrounding water when the ions are located in those regions. In contrast, the probability density curves for water molecules in the bulk-like region show that the local tetrahedral structure of surrounding water is noticeably disrupted by the presence of hydronium, even when the ion itself resides at the interface or within the intermediate region. The corresponding distributions shift toward lower q_tetra values, indicating a transition to a more disordered structure relative to the neutral water system. In comparison, when hydroxide is present, the tetrahedral ordering of nearby water molecules remains largely preserved in the bulk-like region regardless of the ion's location. A slight shift toward lower q_tetra values is observed when hydroxide occupies the bulk-like region, but this effect is less pronounced than for hydronium. These results suggest that hydroxide induces only minor perturbations to the intrinsic layering and hydrogen-bond network, whereas hydronium produces substantial structural disruption extending even into bulk-like water. These results provide a mechanistic explanation for why both ions are found with very low probability in the bulk region.

Fig. 6 Probability density distributions of the tetrahedral order parameter for water molecules in the bulk-like region at a confinement distance of 20 Å as a function of the location of the water ion: (top) hydronium, (bottom) hydroxide. The plots compare three cases: (i) when the ion resides in the interfacial region, (ii) in the intermediate region, and (iii) in the bulk-like region. Each curve is shown relative to the probability density distribution of the tetrahedral order parameter for bulk-like water molecules in the neutral system. For each frame, water molecules were assigned to one of the three regions defined in Fig. S4, and the position of the water ion was evaluated over the full trajectory to determine which region it occupied at each timestep. The tetrahedral order parameter distributions were then computed using all water molecules within the corresponding region across the entire trajectory, and each distribution was normalized so that its total integrated area equals 1.

It is clear that hydronium's tendency to preserve its bulk-like solvation structure under confinement (as indicated by Fig. 2) leads to pronounced disruption of the surrounding water molecules when the ion is located away from the interface. In contrast, when hydronium resides at the interface, favorable surface interactions stabilize it (as indicated by Fig. 5). Although the ion still perturbs the nearby water molecules even at the surface, the energetic benefit of interfacial stabilization outweighs the accompanying structural disruption. Hydroxide behaves inversely because it better accommodates the structured hydrogen-bond network in the confined interior, it can reside in the intermediate and bulk-like regions without significantly altering the surrounding water molecules’ tetrahedral ordering. Thus, the interfacial affinity of hydronium and the partial interior preference of hydroxide arise not solely from their individual solvation structures, but the ion's effect on the surrounding water molecules’ solvation structure under confinement. Hydronium's disruption of interior ordering reinforces its preference for the surface, while hydroxide integrates more compatibly within the ordered water layers.

4 Conclusions

In this work, we investigated the behavior of water ions at the confined hBN–water interface as a function of confinement width. Our results reveal a pronounced preference of hydronium for the interface across all confinement regimes, while hydroxide also exhibits a substantial (but weaker) interfacial affinity. The electronegativity difference between nitrogen and boron in hBN creates a surface dipole. Nitrogen sites can act as hydrogen bond acceptors for hydroxide, partially offsetting the dehydration penalty and enhancing OH⁻ interfacial propensity compared to graphene nanochannels,³⁸ though still less than hydronium. The enhanced interfacial polarization response of hBN therefore modifies the delicate balance of interactions governing hydroxide behavior, enabling partial stabilization near the surface without full disruption of its hydration structure.

Overall, our findings highlight the intricate interplay among hydrogen bonding, orientational ordering, interfacial water structuring, and electrostatic polarization that determines the stability of hydronium and hydroxide at the hBN–water interface. By employing machine learning-based molecular dynamics simulations with first-principles accuracy, we accessed time and length scales far beyond the reach of conventional AIMD, enabling a statistically converged, microscopic description of ion behavior under 2D confinement. The mechanistic insights uncovered here provide a unifying framework for understanding ion-specific interactions at polar 2D interfaces and offer predictive guidance for future experimental and theoretical efforts in nanofluidics, membrane science, and interfacial electrochemistry. This study provides a detailed picture of protonic defects at the hBN–water interface as a foundation for interpreting experimental data and advancing our fundamental understanding of ions at interfaces. Our prediction of slightly higher proton propensity to the hBN–water interface opens up the possibility for technological innovations in nanofluidics, heterogeneous solid–liquid catalysis, and other critical domains that rely on proton-mediated processes.

Author contributions

M. L. B., M. C. A. and T. A. P. designed the research. M. L. B. performed the research. M. L. B. and M. C. A. analyzed the data. M. L. B., C. H. P., M. C. A., and T. A. P. wrote the paper.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data for this article, including input files for CP2K single point and md simulations, input files for training the machine learning potential and input files for running the production MD simulations with ase are included as supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d5fd00126a.

Acknowledgements

This work was supported as part of the Center for Enhanced Nanofluidic Transport (CENT), an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award DE-SC0019112. Work at the Lawrence Livermore National Laboratory was performed under the auspices of the U.S. Department of Energy under Contract DE-AC52-07NA27344. Computing support was provided by the Lawrence Livermore National Laboratory Institutional Computing Grand Challenge program.

References

1. Q. Xie, M. A. Alibakhshi, S. Jiao, Z. Xu, M. Hempel, J. Kong, H. G. Park and C. Duan, Nat. Nanotechnol., 2018, 13, 238–245.
2. Z. Liu, Y. Zhai, K.-G. Zhou and L.-Y. Chu, Adv. Membr., 2021, 1, 100012.
3. S. Wang, L. Yang, G. He, B. Shi, Y. Li, H. Wu, R. Zhang, S. Nunes and Z. Jiang, Chem. Soc. Rev., 2020, 49, 1071–1089.
4. N. R. Aluru, F. Aydin, M. Z. Bazant, D. Blankschtein, A. H. Brozena, J. P. de Souza, M. Elimelech, S. Faucher, J. T. Fourkas and V. B. Koman, et al., Chem. Rev., 2023, 123, 2737–2831.
5. D. Cohen-Tanugi and J. C. Grossman, Nano Lett., 2012, 12, 3602–3608.
6. S. P. Surwade, S. N. Smirnov, I. V. Vlassiouk, R. R. Unocic, G. M. Veith, S. Dai and S. M. Mahurin, Nat. Nanotechnol., 2015, 10, 459–464.
7. H. Li, Z. Song, X. Zhang, Y. Huang, S. Li, Y. Mao, H. J. Ploehn, Y. Bao and M. Yu, Science, 2013, 342, 95–98.
8. S. Joseph, R. J. Mashl, E. Jakobsson and N. R. Aluru, Nano Lett., 2003, 3, 1399–1403.
9. H. G. Park and Y. Jung, Chem. Soc. Rev., 2013, 43, 565–576.
10. Q. Xie, F. Xin, H. G. Park and C. Duan, Nanoscale, 2016, 8, 19527–19535.
11. K. Zhang, Y. Feng, F. Wang, Z. Yang and J. Wang, J. Mater. Chem. C, 2017, 5, 11992–12022.
12. Q. Cai, D. Scullion, W. Gan, A. Falin, S. Zhang, K. Watanabe, T. Taniguchi, Y. Chen, E. J. G. Santos and L. H. Li, Sci. Adv., 2019, 5, eaav0129.
13. L. H. Li, J. Cervenka, K. Watanabe, T. Taniguchi and Y. Chen, ACS Nano, 2014, 8, 1457–1462.
14. M. Carlin, S. Sosa, V. J. González, A. Tubaro, E. Vázquez, M. Prato and M. Pelin, J. Hazard. Mater., 2025, 494, 138449.
15. G. Tocci, L. Joly and A. Michaelides, Nano Lett., 2014, 14, 6872–6877.
16. T. Ohto, H. Tada and Y. Nagata, Phys. Chem. Chem. Phys., 2018, 20, 12979–12985.
17. F. L. Thiemann, C. Schran, P. Rowe, E. A. Müller and A. Michaelides, ACS Nano, 2022, 16, 10775–10782.
18. A. Seal and A. Govind Rajan, Nano Lett., 2021, 21, 8008–8016.
19. A. K. Verma and B. B. Sharma, Surf. Interfaces, 2024, 51, 104646.
20. A. Govind Rajan, M. S. Strano and D. Blankschtein, Nano Lett., 2019, 19, 1539–1551.
21. Z. Li, Y. Li, Y. Yao, F. Aydin, C. Zhan, Y. Chen, M. Elimelech, T. A. Pham and A. Noy, ACS Nano, 2020, 14, 6269–6275.
22. C. Zhan, M. R. Cerón, S. A. Hawks, M. Otani, B. C. Wood, T. A. Pham, M. Stadermann and P. G. Campbell, Nat. Commun., 2019, 10, 4858.
23. Z. Li, A. T. Hall, Y. Wang, Y. Li, D. O. Byrne, L. R. Scammell, R. R. Whitney, F. I. Allen, J. Cumings and A. Noy, Sci. Adv., 2024, 10, eado8081.
24. N. Kavokine, M.-L. Bocquet and L. Bocquet, Nature, 2022, 602, 84–90.
25. Y. Yu, J. Fan, A. Esfandiar, Y. Zhu, H. Wu and F. Wang, J. Phys. Chem. C, 2019, 123, 1462–1469.
26. Q. Yang, P. Z. Sun, L. Fumagalli, Y. V. Stebunov, S. J. Haigh, Z. W. Zhou, I. V. Grigorieva, F. C. Wang and A. K. Geim, Nature, 2020, 588, 250–253.
27. Y. Xue, Y. Xia, S. Yang, Y. Alsaid, K. Y. Fong, Y. Wang and X. Zhang, Science, 2021, 372, 501–503.
28. S. Goutham, A. Keerthi, A. Ismail, A. Bhardwaj, H. Jalali, Y. You, Y. Li, N. Hassani, H. Peng, M. V. S. Martins, F. Wang, M. Neek-Amal and B. Radha, Nat. Nanotechnol., 2023, 18, 596–601.
29. Y. Yu, J. Fan, J. Xia, Y. Zhu, H. Wu and F. Wang, Nanoscale, 2019, 11, 8449–8457.
30. F. Aydin, A. Moradzadeh, C. L. Bilodeau, E. Y. Lau, E. Schwegler, N. R. Aluru and T. A. Pham, J. Chem. Theory Comput., 2021, 17, 1596–1605.
31. T. A. Pham, S. M. G. Mortuza, B. C. Wood, E. Y. Lau, T. Ogitsu, S. F. Buchsbaum, Z. S. Siwy, F. Fornasiero and E. Schwegler, J. Phys. Chem. C, 2016, 120, 7332–7338.
32. L. Fumagalli, A. Esfandiar, R. Fabregas, S. Hu, P. Ares, A. Janardanan, Q. Yang, B. Radha, T. Taniguchi and A. K. Geim, et al., Science, 2018, 360, 1339–1342.
33. R. Wang, M. Souilamas, A. Esfandiar, R. Fabregas, S. Benaglia, H. Nevison-Andrews, Q. Yang, J. Normansell, P. Ares, G. Ferrari, A. Principi, A. K. Geim and L. Fumagalli, Nature, 2025, 646, 606–610.
34. Y. Wang, H. Luo, X. R. Advincula, Z. Zhao, A. Esfandiar, D. Wu, K. D. Fong, L. Gao, A. S. Hazrah, T. Taniguchi, C. Schran, Y. Nagata, L. Bocquet, M.-L. Bocquet, Y. Jiang, A. Michaelides and M. Bonn, J. Am. Chem. Soc., 2025, 147, 30107–30116.
35. B. Grosjean, M.-L. Bocquet and R. Vuilleumier, Nat. Commun., 2019, 10, 1656.
36. B. Grosjean, C. Pean, A. Siria, L. Bocquet, R. Vuilleumier and M.-L. Bocquet, J. Phys. Chem. Lett., 2016, 7, 4695–4700.
37. A. Hosseini, A. Masoud Yarahmadi, S. Azizi, A. H. Korayem and R. Savary, Phys. Chem. Chem. Phys., 2024, 26, 21841–21849.
38. X. R. Advincula, K. D. Fong, A. Michaelides and C. Schran, ACS Nano, 2025, 19, 17728–17737.
39. L. Scalfi, L. Lehmann, A. P. dos Santos, M. R. Becker and R. R. Netz, J. Chem. Phys., 2024, 161, 144701.
40. Y.-L. S. Tse, C. Chen, G. E. Lindberg, R. Kumar and G. A. Voth, J. Am. Chem. Soc., 2015, 137, 12610–12616.
41. P. Zhang, M. Feng and X. Xu, ACS Phys. Chem. Au, 2024, 4, 336–346.
42. S. Luo, R. P. Misra, A. Sam, T. Frömbgen, B. Kirchner and D. Blankschtein, ACS Nano, 2025, 19, 27500–27514.
43. G. He, Q. Pan, Z. Yuan, L. F. Villalobos, J.-H. Park, H.-Y. Chi, M. Kuehne, Y. Zeng, Y.-M. Tu, J. Zhao, J. Kong, Z. Jiang, K. V. Agrawal, D. Blankschtein and M. S. Strano, Angew. Chem., Int. Ed., 2025, 64, e202514288.
44. S. Luo, R. P. Misra and D. Blankschtein, ACS Nano, 2024, 18, 1629–1646.
45. T. Dufils, C. Schran, J. Chen, A. K. Geim, L. Fumagalli and A. Michaelides, Chem. Sci., 2024, 15, 516–527.
46. I. Batatia, D. P. Kovács, G. N. C. Simm, C. Ortner and G. Csányi, MACE: Higher Order Equivariant Message Passing Neural Networks for Fast and Accurate Force Fields, arXiv, 2023, preprint, arXiv:2206.07697 [stat],  DOI:10.48550/arXiv.2206.07697.
47. T. D. Kuhne, M. Iannuzzi and M. Del Ben, et. al., J. Chem. Phys., 2020, 152, 194103.
48. J. P. Perdew, Phys. Rev. Lett., 1996, 77, 3865–3868.
49. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
50. B. Cheng, E. A. Engel, J. Behler, C. Dellago and M. Ceriotti, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 1110–1115.
51. J. VandeVondele and J. Hutter, J. Chem. Phys., 2007, 127, 114105.
52. I. Batatia, P. Benner, Y. Chiang, A. M. Elena, D. P. Kovács, J. Riebesell, X. R. Advincula, M. Asta, M. Avaylon, W. J. Baldwin, F. Berger, N. Bernstein, A. Bhowmik, F. Bigi, S. M. Blau, V. Cărare, M. Ceriotti, S. Chong, J. P. Darby, S. De, F. Della Pia, V. L. Deringer, R. Elijošius, Z. El-Machachi, E. Fako, F. Falcioni, A. C. Ferrari, J. L. A. Gardner, M. J. Gawkowski, A. Genreith-Schriever, J. George, R. E. A. Goodall, J. Grandel, C. P. Grey, P. Grigorev, S. Han, W. Handley, H. H. Heenen, K. Hermansson, C. H. Ho, S. Hofmann, C. Holm, J. Jaafar, K. S. Jakob, H. Jung, V. Kapil, A. D. Kaplan, N. Karimitari, J. R. Kermode, P. Kourtis, N. Kroupa, J. Kullgren, M. C. Kuner, D. Kuryla, G. Liepuoniute, C. Lin, J. T. Margraf, I.-B. Magdău, A. Michaelides, J. H. Moore, A. A. Naik, S. P. Niblett, S. W. Norwood, N. O'Neill, C. Ortner, K. A. Persson, K. Reuter, A. S. Rosen, L. A. M. Rosset, L. L. Schaaf, C. Schran, B. X. Shi, E. Sivonxay, T. K. Stenczel, C. Sutton, V. Svahn, T. D. Swinburne, J. Tilly, C. van der Oord, S. Vargas, E. Varga-Umbrich, T. Vegge, M. Vondrák, Y. Wang, W. C. Witt, T. Wolf, F. Zills and G. Csányi, J. Chem. Phys., 2025, 163, 184110.
53. G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126, 014101.
54. A. Hjorth Larsen, J. Jørgen Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. Bjerre Jensen, J. Kermode, J. R. Kitchin, E. Leonhard Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. Bergmann Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiøtz, O. Schütt, M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. Zeng and K. W. Jacobsen, J. Phys.: Condens. Matter, 2017, 29, 273002.
55. G. J. Martyna, M. E. Tuckerman, D. J. Tobias and M. L. Klein, Mol. Phys., 1996, 87, 1117–1157.
56. S. Shekarforoush, H. Jalali, M. Yagmurcukardes, M. V. Milošević and M. Neek-Amal, Phys. Rev. B, 2020, 102, 235406.
57. F. L. Hirshfeld, Theor. Chim. Acta, 1977, 44, 129–138.
58. A. Luzar and D. Chandler, Nature, 1996, 379, 55–57.
59. E. Duboue-Dijon and D. Laage, J. Phys. Chem. B, 2015, 119, 8406–8418.

---