Molecular picture of electric double layers with weakly adsorbed water

Water molecules adsorb on electrode surfaces when the solid and liquid phases come into contact with each other.¹ It has been shown that the adsorbed water molecules play an important role in many sustainable chemical technologies, such as (photo)electrocatalytic water splitting,² water desalination,³ and water harvesting.⁴ The water adsorption energy E_ads is a key physical quantity that can reflect the interaction between electrolyte and electrode materials.⁵ Therefore, unraveling the relationship between E_ads and the surface properties of electrode materials is of paramount importance.

When in contact with liquid electrolyte under working conditions, the electrode materials will attract counter-ions from the solution and form the electric double layers (EDLs) that strongly modulate the interfacial reactivity and wetting. EDLs can be developed either at metal electrodes beyond the potential of zero charge⁶ or at metal-oxide electrodes outside the point of zero charge. In the latter case, the metal-oxide surfaces are positively or negatively electrified by the adsorption of protons or hydroxyl groups, which depend on the pH of the environment.⁷ 

The well-known classical Gouy–Chapman–Stern (GCS) model has been used to describe the structure of EDLs at oxide surfaces, in which EDLs consist of a space charge layer, a Helmholtz layer, and a Gouy–Chapman layer.^1,8 The Helmholtz layer is formed by the interfacial proton charges and counter-ions from the electrolyte, which spans for 3–5 Å. The microscopic structure, in particular the fully solvated surface and the interaction and orientation of adsorbed molecular or dissociated water molecules, determines the electrochemical performances, i.e., the thermodynamics, kinetics, and product selectivity of electrode materials. Therefore, the description of the relationship between water adsorption energy E_ads and EDLs is crucial.^9–13 

Thanks to the development of advanced in situ experimental technologies and theoretical studies, it was not until recently that the structure and properties of interfacial EDLs started to be revealed. For example, electrochemical x-ray photoelectron spectroscopy (EC-XPS) applies to the characterization of the surface structure of the electrode,¹⁴ sum-frequency vibrational spectroscopy retrieves the vibrational spectra of the bonded water layer and the ion diffuse layer,^15, in situ Raman spectroscopy investigates the interfacial water on single-crystal surfaces,^16,17 the dual scale atomic force microscopy (AFM) spectroscopy technique provides direct access to the facet-dependent local surface charge density and the potential drop across the diffuse part of the EDL,¹⁸ and operando scanning electrochemical probe microscopy can correlate electrochemical activity with changes in surface properties.¹⁹ Moreover, various computational methods, e.g., quantum mechanics/implicit solvation model (QM/ISM),^20,21 quantum mechanics/molecular mechanical (QM/MM),^22,23 density functional theory-based molecular dynamics simulation (DFTMD),^24–27 and machine learning-accelerated MD simulations,^28–30 have been employed to simulate the electrochemical interfaces and describe the microscopic structure and the dynamics of water and ions near the interfaces.

Despite the tremendous experimental and theoretical studies that have been devoted to this topic, a molecular-level understanding of the effect of water adsorption energy E_ads on microscopic structure and Helmholtz capacitance within EDLs remains elusive, owing to its complexity and difficulty in being probed with sufficient resolution.^31–33 Recently, by combining the modern theory of polarization and the constant electric displacement D Hamiltonian developed by Stengel, Spaldin, and Vanderbilt,³⁴ the finite-field DFTMD method has become an attractive technique for modeling EDLs at electrified solid–liquid interfaces.^35–38 With the same amount but opposite types of proton charges, the two slabs of the simulated oxide surface are charged at a fixed chemical composition of the supercell. In the finite-field DFTMD simulations, the average Helmholtz capacitance C_H is obtained based on the supercell polarization P_z.³⁹ Note that the ensemble averaged supercell polarization $Pz$ converges rapidly by the application of constant electric displacement simulations, which makes the computation of C_H accessible within the timescale of DFTMD simulation.³⁶ 

In the previous work on applying finite-field DFTMD for modeling EDLs at both rutile TiO₂(110)/NaCl³⁶ and SnO₂(110)/NaCl³⁷ interfaces, it was found that water chemisorption is a key factor in accounting for the asymmetric EDLs seen at aqueous interfaces of oxide. To further quantify the relationships between E_ads, molecular structure, and Helmholtz capacitance of EDLs on metal oxides, one of the most promising n-type oxide bismuth vanadates, BiVO₄, is selected here as an example of an oxide electrode with weakly chemisorbed water. BiVO₄ has been applied to photoelectrochemical water splitting and has garnered considerable attention due to its low-cost, stability, nontoxicity, ideal bandgap (∼2.4 eV), and high electron–hole separation efficiencies (over 70%).^40,41 In addition, the E_ads of BiVO₄ was reported to be −0.5 eV,⁴² which is much smaller than that of TiO₂(110) (−0.87 eV⁴³) and SnO₂(110) (−1.5 eV³⁷). Moreover, it can be seen as a system with weakly chemisorbed water.

In the following, we first validate the DFTMD simulations of BiVO₄(010)/NaCl electrolyte interfaces by computing the corresponding surface pK_as and comparing them to experiments and the results of TiO₂(110) and SnO₂(110) surfaces. It is followed by the structural and dynamical analyses of water orientation, the local coordination environment, and counter-ions in the protonic double layers of BiVO₄(010). Finally, we link the degree of asymmetry in the computed Helmholtz capacitance and the adsorption energy of water, which shows a surprising positive correlation at both metal oxide and graphene-like surfaces.

The acidity constants and EDLs at BiVO₄(010)/NaCl electrolyte interfaces are calculated at a fixed chemical composition of the supercell. The BiVO₄(010) surface was modeled by symmetric periodic slabs of eight-layers with lateral dimensions of a 2 × 2 surface cell. Note that the slab size was tested by the convergence calculations of E_ads and bandgap (see Table S1 of the supplementary material). The two symmetric slabs were separated by an electrolyte space of 30.5 Å, giving an orthorhombic supercell of dimensions 10.39 × 10.18 × 53.25 Å³. A molality of 2.7 mol/kg (the concentration of the electrolyte used in industry) for NaCl electrolyte was simulated in our systems. The electrolyte space between the BiVO₄(010) slabs was filled with 102 water molecules, 5 Na⁺ and 5 Cl⁻ in each unit cell.

All the DFTMD simulations are carried out using the freely available program package CP2K/Quickstep^44,45 with the PBE functional.⁴⁶ Goedecker–Teter–Hutter (GTH) pseudopotentials are employed to represent the core electrons.^47,48 The atomic basis sets for the valence electrons were the standard short-ranged double-ζ basis functions, with one set of polarization functions (DZVP) optimized for molecular systems.⁴⁹ The plane wave cutoff for electron density expansion is set at 400 Ry. The target accuracy for the SCF convergence was 3 × 10⁻⁶ a.u. The NVT ensemble is used for MD propagation with a time step of 0.5 fs, and the equilibrium temperature for all the simulations is kept to 330 K by using the Nose–Hoover thermostat.

For further details on acidity calculations and finite-field DFTMD simulations of BiVO₄(010)-NaCl interfaces, see the supplementary material.

The classical MD simulations of the graphene/NaCl interfaces were conducted using the LAMMPS software package.⁵⁰ The electrode material is composed of two graphene layers with a surface area of 51.1 × 49.2 Å². The two electrodes were separated by a vacuum region of 13.2 Å, giving the total length of the simulated box to be 70 Å. The NaCl electrolyte is comprised of 4100 water molecules and 74 pairs of NaCl ions, giving a concentration of 1 M. The computational details of classical MD simulations are available in the supplementary material.

When the BiVO₄ electrode is immersed in an aqueous electrolyte, the protonation or deprotonation of surface groups depends on their pK_a values and pH in the solution. There is a special pH at which the concentration of H⁺ on the surface equals that of OH⁻, i.e., the point of zero charge (PZC).⁶² PZC can be measured using the potentiometric titration method and electrochemical in situ spectroscopy technologies.^53,63 However, the intrinsic pK_a values of individual surface groups are typically unable to be characterized by experiments. Herein, we first calculate the acidity constants of the individual adsorption sites on the BiVO_a(010)/NaCl interface from DFTMD simulations and free energy perturbation theory by applying harmonic restraining potentials V_rH to all O–H bonds (see Table S2 of the supplementary material).^64,65

In this work, we have used the full reaction scheme in the calculation of the surface acidity constants developed previously,^65,67,68 as illustrated in Fig. 1. The deprotonation of the Bi_7fOH₂ group was coupled to a simultaneous protonation of H₂O in the bulk electrolyte of the same inhomogeneous model system [Fig. 1(a), Reaction S2 of the supplementary material]. Similarly, the deprotonation process of a BiVO_brH⁺ group was accompanied by another protonation of H₂O in the bulk electrolyte [Fig. 1(b), Reaction S3 of the supplementary material]. The reverse of this reaction indicates a bridging oxygen site receiving a proton from the acidic bulk electrolyte. The corresponding deprotonation free energies ΔA₁ and ΔA₂ were calculated from a coupling parameter integral of the corresponding vertical energy gap (see Figs. S1 and S2 of the supplementary material). The pK_a values we found for the two active surface hydroxyl groups are 6.2 and 0.8, respectively, giving the PZC value of 3.5 (Table I), in excellent agreement with the experimental characterization 2.5–4.5^51–53 and previous theoretical works.⁶⁶ 

Since the H-bond network has a significant effect on pK_a values, both the coordination numbers and the length of the H-bond can lead to a vital variation in the picture of surface acidity. To further interpret the effect of water adsorption energy E_ads on pK_a shifts, the RDFs between oxygen atoms of surface groups and hydrogen atoms in electrolyte are shown in Fig. 1. The comparison of H-bond coordination numbers of BiVO₄ with those of TiO₂(110) and SnO₂(110) surfaces is shown in Table II. Note that the H-bond coordination numbers are obtained by integrating the first peaks of RDFs in Figs. 1(c) and 1(d), and the length of the H-bond for oxygen atoms on the BiVO₄ surface is determined by the corresponding position of the first peak. The E_ads of BiVO₄(010) surface was calculated to −0.58 eV by one monolayer of adsorbed molecular water (Table I).

It is known that the surface metal–O bond strength is manifested by the corresponding bond-length.⁶⁹ Because of the weakest interaction between the adsorbed water molecules and BiVO₄ surface sites, it is found that the bond length of Bi_7f and oxygen in the adsorbed water molecules is 2.6 Å, which is significantly larger than that of TiO₂(110) and SnO₂(110) surfaces (∼2.0 Å).⁴³ On the other hand, the interaction between adsorbed water and solvent water in an electrolyte increases when the strength of water adsorption weakens. It is indeed found to have a higher H-bond coordination number of 0.97 for the Bi_7fOH₂ group than that of the TiOH₂ (0.6) and SnOH₂ (0.3) groups (Table II). In the case of the deprotonated group, it is found that the length of the H-bond of Bi_7fOH⁻ [first peak of the blue line in Fig. 1(c)] is significantly smaller than that of Bi_7fOH₂ [first peak of the blue dashed line in Fig. 1(c)], while the length of the H-bond of TiOH⁻ is the same as that of TiOH₂.⁴³ One would expect that a shorter H-bond (stronger H-bonding) formed by Bi_7fOH⁻ makes the deprotonated structure more stable, and more H-bonds (stronger H-bonding) formed on the Bi_7fOH₂ site make the protonated state more stable. Indeed, these two effects appear to cancel out each other and leave the resulting pK_a1 of the BiVO₄ surface similar to that of TiO₂ (Table I). Meanwhile, it was found that BiVO_brH⁺ has a significantly lower coordination number 0.0 than that of the Ti₂O_brH⁺ (1.0) and Sn₂O_brH⁺ (1.0) groups (Table II). This difference is much bigger than the lowering of the coordination number for deprotonated state BiVO_br in the oxide series. Therefore, the protonated state of BiVO_brH⁺ is much less stable, and the resulting pK_a2 of the BiVO₄ surface is the smallest.

Based on the values of pK_a, we obtained the dissociation free energy of adsorbed water ΔA_diss on the BiVO₄(010) surface according to Eq. (6), giving a value of 0.32 eV. There have always been different views on the dissociation of chemisorbed water at the BiVO₄ interface. Previous theoretical studies had not reported or observed the dissociation of adsorbed water on the BiVO₄(010) surface until Wang and Galli found that the dissociation of water does occur on the BiVO₄(010) surface using ambient pressure resonant photoemission spectroscopy (AP-ResPES) measurements and first-principles simulations.² It is worth noting that water dissociation on the pristine rutile TiO₂(110) surface has been controversial for both the experiment and theory, with the dissociation free energy of adsorbed water ranging from 0.23 eV⁵⁴ to 0.45 eV⁷⁰ with density functional theory (DFT) studies. Compared to TiO₂, BiVO₄ has a slightly higher ΔA_diss (see Table I). Therefore, it is plausible for the spontaneous dissociation of adsorbed water on the BiVO₄ surface with a sufficient timescale, although we did not observe it during the course of our DFTMD simulations at PZC.

Among rutile TiO₂(110) and SnO₂(110) and BiVO₄(010) surfaces, the SnO₂(110) surface has the strongest water adsorption energy of −1.5 eV with DFT studies³⁷ (see Table I). Therefore, one would expect water molecules to have a stronger tendency for dissociative adsorption and ΔA_diss to be smaller. This is indeed the case, as shown in Table I.

In previous reports, it has been known that the local surface potentials and photoelectrochemical performance depend strongly on the pH of the ambient electrolyte^40,53 and the interaction of H₂O with the BiVO₄ surface.^2,71 To further understand how the chemisorbed water response to surface potential affects the molecular structure and capacitance of EDLs, we next simulate the electrified BiVO₄/NaCl interface with finite-field DFTMD simulations (Fig. S3 of the supplementary material).

At low pH, the BiVO₄ surface is positively charged by adsorbing protons, and the surface charge density σ = 30.3 μC/cm². At high pH, the BiVO₄ surface is negatively charged by the desorption proton from terminal Bi_7fOH₂, and σ = −30.3 μC/cm². Detailed descriptions of the computational setup are given in the supplementary material. When looking at the structure of EDLs, it is found that the electrified BiVO₄/NaCl interfaces are very different from those of rutile oxide TiO₂ and SnO₂, due to weaker water adsorption energy.

At BiVO₄(010)/NaCl interfaces, the dynamics of adsorption and desorption processes can be presented by the time evolution of the Bi–O_w distance between surface Bi and oxygen in water, as shown in Figs. 3(a)–3(c). When the Bi–O_w distance goes from ∼5 to ∼2.6 Å, it indicates the adsorption process of second layer water. When the Bi–O_w distance goes from ∼2.6 to ∼5 Å, it indicates the desorption process of the first layer of water. Interestingly, it is found that the dynamics of adsorption and desorption processes have a big difference between the first layer of water and the second layer of water at different charged BiVO₄(010)/NaCl interfaces.

At PZC, the adsorption process [black dashed line in Fig. 3(a)] is accompanied by the desorption process [gray dashed line in Fig. 3(a)] between the second layer of water and the first layer of water. Furthermore, it is observed that more than one water molecule adsorbed to the same Bi_7f site at the positively charged BiVO₄/NaCl interface [Fig. 3(b)]. That means the adsorption of first layer water is enhanced by the electric field in EDL at low pH. However, the adsorbed first layer of water partly desorbed into the near-surface at the negatively charged BiVO₄/NaCl interface [Fig. 3(c)]. The results indicate that the adsorption of first layer water is weakened by the electric field in EDL at high pH.

The RDFs between Bi_7f sites and oxygen atoms of electrolyte water at different charged BiVO4(010)/NaCl interfaces are described in Fig. 3(d). The average number of adsorbed water molecules for each Bi_7f site is obtained by integrating the first peak of the Bi_7f–O RDFs, giving 1.0, 1.3, and 0.6 at PZC, respectively, at positively charged and negatively charged interfaces. That means the amount of adsorbed water at the BiVO₄ surface depends on its surface charge density, while it remains unchanged at the TiO₂ and SnO₂ surfaces. Although the adsorption and desorption processes can happen dynamically on the negatively charged side of TiO₂, the number of adsorbed water is the same as in the cases of PZC and the positively charged side.³⁶ These results give us a molecular-level understanding of the structure and dynamics of metal oxide electrodes with weakly chemisorbed water, which is different from the oxides with strong water adsorption energy.

It has been shown that the mechanism of the water oxidation reaction on the BiVO₄ surface depends on the pH of the electrolyte.⁷² At low pH, the water oxidation reaction is initiated mainly through the dehydrogenation of adsorbed water molecules, which involves the H₂O^·+ radical cation as an intermediate.⁷³ At high pH, the reaction is initiated through the oxidation of OH⁻ (OH⁻/OH^·).⁷⁴ Based on observations here, the Bi_7f–(OH₂)₂ structure may be more favorable to the formation of H₂O^·+ in acidic conditions. However, the desorption of terminally adsorbed water is a disadvantage to generate H₂O^·+ at alkaline conditions. This leaves the opportunity for Bi_7fOH⁻ to combine hole and formate OH^· at high pH. This is in line with the role of negative surface charge in hole selectivity reported by recent experiments^75,76 and the ionized oxygen vacancy proposed by theoretical calculations.⁷⁷ 

In addition, such a weak interaction between interfacial water and BiVO₄(010) surface groups gives the compensated counterions in the electrolyte the chance to have a direct and strong interaction with surface sites. Therefore, the inner sphere coordination formed at the BiVO₄(010)/NaCl surface [see Figs. 4(a) and 4(b)]. The width of the Helmholtz layer d_EDL is estimated at about 2.4 Å at low pH and 2.7 Å at high pH, respectively [Fig. 4(c)]. Note that the values of d_EDL were obtained from the surface normal projected distance between the counterions and the nuclei of O_br sites. In addition, it is found that the adsorbed Cl⁻ on the Bi_7f site at a positively charged surface [Fig. 4(d)] makes the width of d_EDL smaller than that of the negatively charged surface. In the reported synchronous illumination XPS (SI-XPS) study, the Cl⁻ concentrates and bonds to Bi on the surface, which can facilitate H₂O molecule activation as well as promote carrier charge separation and accelerate the water oxidation reaction.^78,79 This may be explained by the adsorption of Cl⁻, which shortens the width of d_EDL and has a positive effect on Helmholtz’s capacitance, as seen here.

When looking at Table III and Fig. 2 together, one can see the shift in the orientation peak of water at positively charged surface Δθ⁺ is the largest in the case of BiVO₄, which correlates well with $CH+$ (BiVO₄) being the highest in the series. On the other hand, the shift in the orientation peak of water at negatively charged surface Δθ⁻ is the smallest in the case of BiVO₄, which also correlates with the lowest capacitance value $CH−$(BiVO₄) among these three oxide systems. However, this trend is not clear when just comparing TiO₂ and SnO₂. This means this qualitative picture based on water orientation alone is not sufficient, and the incorporation of the different water adsorption strengths and their effects on polarization is needed in order to be quantitative.

In previous work, it has been shown that water chemisorption would explain asymmetric EDLs seen at aqueous interfaces of oxide.^36,37 In this work, it is found that EDLs are indeed nearly symmetric on the BiVO₄(010) surface, which possesses the weakest chemisorption of water molecules (Table III). When looking at the capacitance asymmetry of EDLs in metal oxide systems, a conspicuous monotonic increase is found in ΔC_H when increasing the water adsorption energy of E_ads, as shown in Fig. 5(c).

To verify this relationship between the capacitance asymmetry ΔC_H and E_ads at other types of aqueous interfaces, i.e., the metallic electrode with weakly physisorbed water, the graphene/NaCl electrolyte interface was selected and simulated with classical force fields and constant-potential MD (see the supplementary material and Fig. S6).⁸⁰ It is worth noting that the nature of EDL at the metallic interface is quite different from that of the (insulating/semiconducting) oxide interface; the former is formed by the electronic charge [Fig. 5(b)], while the latter is formed by the proton charge [Fig. 5(a)].

The adsorption energy was calculated from a single water molecule with a two-leg structure adsorbed on the surface of a two-layer graphene electrode (see Fig. S7 of the supplementary material). The strength of the water adsorption energy E_ads of graphene electrodes was changed by adjusting the energy parameter ɛ of the Lennard-Jones potential between water molecules and graphene. The E_ads in the classical MD simulations with energy parameters of 0.3ɛ, 0.7ɛ, ɛ, 1.5ɛ, and 2.0ɛ were calculated as 0.05, 0.08, 0.10, 0.14, and 0.17 eV, respectively. In the classical MD simulations, the Helmholtz capacitance (C_H) is calculated from the definition of C_H = dσ/d(Δϕ), and Δϕ is the electrostatic potential difference between the top graphene layer and bulk electrolyte, as shown in Fig. S8 of the supplementary material. The calculated $CH+$ and $CH−$ for electrified graphene/NaCl electrolyte interfaces with different adsorption capabilities are described in Fig. S9 of the supplementary material.

Based on the classical MD simulations, it is found that the ΔC_H shows a similar monotonic increment as a function of E_ads at electrified graphene/NaCl electrolyte interfaces [Fig. 5(d)]. This result indicates that the positive correlation that we found on oxide surfaces is also seen in graphene-like systems.

The values of the asymmetry factor γ are shown on the right side of Figs. 5(c) and 5(d). Despite the fact that the absolute values of C_H are dramatically different between metal-oxide and graphene-like surfaces, the asymmetry factor γ shows similar values. This indicates that the relationship between E_ads and γ can be rather general across different types of electrode materials, although different scaling factors apply depending on the nature of EDLs. Furthermore, it is worth noting that the water adsorption energy at the graphene surface, as determined from highly accurate electronic structure calculations, is about 0.1 eV,⁸¹ which corresponds to a high value in the asymmetry factor γ of 0.25 in our quantification. This can be the physical reason behind the observed asymmetric response of the contact angle in the electrowetting experiment recently performed on this type of material.⁸² 

In conclusion, by applying finite-field DFTMD simulations, the effect of water adsorption energy E_ads on surface acidity and Helmholtz capacitance has been studied by selecting BiVO₄ as an example of metal oxide electrodes with weakly chemisorbed water. Based on an oxide series composed of BiVO₄, TiO₂, and SnO₂, we investigate the effect of E_ads on the coordination number, the length of the H-bond, and the orientation of chemisorbed water. In particular, it is found that a positive correlation exists between the degree of asymmetry in the Helmholtz capacitance and the strength of the E_ads. This correlation is further verified and extended to graphene-like systems with physisorbed water.

Overall, these observations in this work provide a general molecular picture of an electric double layer with different E_ads. Our findings likely lay the foundation for future investigations of water reactivity and wetting in electrified semiconducting metal oxides and carbonaceous materials.

See the supplementary material for the following additional data: the convergence tests for the slab size, harmonic restraining potentials used in the acidity calculations, calculations of the deprotonation free energy ΔA, finite-field DFTMD simulations of electrified BiVO4(010)–NaCl interfaces, and classical MD simulations of electrified graphene–NaCl electrolyte interfaces.

This project received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation Program (Grant Agreement No. 949012). Mei Jia greatly appreciates the financial support from the Natural Science Foundation of Henan Province (Grant No. 242300420420570), the Key Scientific Research Projects of Colleges and Universities in Henan Province (Grant No. 24A150031), and the International Scientific and Technological Cooperation Projects in Henan Province (Grant No. 232102520020). We thank S. Han for the discussions and inputs during the initiation of this project.

The authors have no conflicts to disclose.

Mei Jia: Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Writing – original draft (equal); Writing – review & editing (equal). Junyi Wang: Data curation (supporting); Formal analysis (supporting); Investigation (supporting); Visualization (supporting). Qixiang Liu: Formal analysis (supporting); Writing – review & editing (supporting). Xiaohui Yang: Data curation (supporting); Formal analysis (supporting); Investigation (supporting); Writing – review & editing (supporting). Chao Zhang: Conceptualization (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding authors upon reasonable request.