Toward rational understanding of the hydrogen evolution polarization curves through multiscale simulations

Domain: physics
Published: 2025-12-20
Editorial outcome: major_revision

Abstract
To properly understand complex heterogeneous electrocatalytic processes, it is essential to consider not only the microscopic electrochemical reaction mechanisms at electrodes but also the macroscopic mass transport effects, which are usually overlooked in current computational electrochemistry. Here, we develop a multiscale simulation framework that integrates grand canonical density functional theory calculations, microkinetic modeling, and a continuum transport model to successfully reproduce the experimental polarization curves of hydrogen evolution reaction on Au(111) across the entire pH range. Our simulations reveal how the potential-dependent evolution of local environments governs reaction mechanisms and behavior, identifying a dynamic balance between local pH increases and decreases in reaction barriers that controls characteristic current plateaus. By classifying polarization curves into five distinct regimes based on their governing factors, we provide comprehensive mechanistic insights into how mass transport influences polarization behavior. Our multiscale simulation framework not only adequately reproduces existing experimental results but also provides a basis for rational optimization of the microenvironment to achieve efficient electrochemical processes.

Introduction
Electrocatalysis has played a vital role in harnessing sustainable energy into solar fuels, where designing high-performance and low-cost electrocatalysts has long been a challenging yet appealing task1–3. Achieving rational catalyst design requires fundamental understanding of electrocatalytic mechanisms, which are typically governed by both intrinsic properties of catalysts and local reaction environments4–7. Currently, the state-of-the-art characterization techniques and spectroscopic measurements greatly benefit the mechanistic studies in electrocatalysis8–16, while the theoretical simulations could further push the understanding to a more microscopic boundary17–26. To understand better or even predict new catalysts for electrocatalysis, several advanced computational methods and models have been developed to accurately calculate the reaction energy profiles for different heterogeneous electrochemical processes, which incorporate influences of electrode potential, pH, solvation, and reaction kinetics27. Based on their fundamental importance, failing to adequately integrate mass transport effects into reaction mechanism investigation may create a significant gap between computational models and realistic electrocatalytic systems.

Specifically, a practical electrochemical reaction will not only involve the microscopic process at the electrode interface, i.e., the adsorption/desorption of reactants/products, bond breaking/formation, and proton-coupled electron transfer (PCET), but also the macroscopic mass transport phenomenon emerging from the collective behavior of numerous atoms and molecules, i.e., diffusion, migration, and convection. In principle, incorporating both microscopic and macroscopic factors into computational electrocatalysis is necessary to achieve a comprehensive understanding of experimental phenomena, but it is also significantly challenging. Currently, although continuum transport models have been employed to simulate mass transport phenomena and reaction environments in macroscopic electrochemical cells28–30, the unquantifiable uncertainties of reaction kinetics associated with traditional Tafel equation-based models still need to be properly addressed31,32. Thus, multiscale simulation is still one of the most promising methods to integrate microscale reaction mechanisms with macroscale transport effects and eventually achieve a complete modeling of the electrochemical cell4,7,33–36.

Indeed, Bell and coworkers reported a pioneering multiscale model that elegantly integrates microscopic reaction mechanisms with macroscopic mass transport phenomena in a device for CO2 electroreduction, thereby explaining how mass transport influences the atomistic reaction mechanisms37. Thereafter, Chan and coworkers achieved qualitative agreement with experimental results because only the diffusion effect was included in their mass transport model38. Then, Liu and coworkers incorporated both diffusion and convection effects into their transient kinetics simulation, enabling an improved reproduction of the experimental polarization behavior of hydrogen evolution reaction (HER) on the Pt(111) surface39. Wang and coworkers also reasonably predicted polarization curves of HER on the Pt(111) surface by building a bilayer explicit water model to mimic the proton transfer from bulk water to the interfacial water40. Close inspections revealed that each work has diverse implementations and varying understandings of mass transport factors in modeling the electrochemical process, making it challenging to conduct a fair comparison and extend the methods conveniently to other new systems. Thus, a simulation workflow with the merits of facile implementation, extensibility, and strict considerations of diffusion, migration, and convection in mass transport is still missing but highly desirable.

In this work, using HER on Au catalysts as the model reaction, we establish a multiscale simulation framework that couples grand canonical density functional theory (GC-DFT), microkinetic modeling (MKM), and continuum transport models, thereby illustrating the relationship between reaction mechanisms and mass transport processes at both microscopic and macroscopic scales. In this framework, the potential-dependent reaction free energies and activation energy barriers for each elementary step were calculated using the GC-DFT method. The reaction rates obtained from microkinetic modeling were then implemented as boundary conditions in the continuum transport model to simulate electrolyte mass transport and determine pH and H2 concentration profiles. Our theoretically simulated results on Au(111) and Pt(111) surfaces demonstrate good agreement with available experimental polarization curves of HER41, enabling us to unravel the origin of the characteristic current plateau. By classifying the polarization curve into five distinct regimes, we elucidate the intricate interplay between reaction kinetics and mass transport.

Results
Multiscale simulation framework
As illustrated in Fig. 1a, the HER on electrodes proceeds through different mechanisms depending on the pH, where the proton sources are typically hydronium ions (H3O+) in acidic HER but water in the alkaline medium. The continuous electrochemical reactions lead to reactant depletion and product accumulation at the electrode, resulting in deviations between the electrode’s local environment and the bulk solution. To quantitatively describe this phenomenon, we integrated the contributions of mass transport with first-principles calculations, thereby establishing a continuum transport model to simulate the relationship between the local environment and reaction rate of HER at the electrode.

Figure 1b and Supplementary Figs. 1-2 illustrate the specific workflow and key components of our multiscale modeling approach. First, the potential-dependent reaction free energies and energy barriers for all elementary steps are calculated using the GC-DFT method, which are then incorporated into the Microkinetic Model (MKM) to obtain the equilibrium constants, the rate constants, surface coverages, and reaction rates (with the bulk pH and H2 concentration as inputs). At this stage, the current density, excluding mass transport effects, can be explicitly calculated. Then, the obtained current density will serve as input for the continuum transport model that resolves the local pH and H2 concentration profiles at the cathode interface. The specific definitions of bulk/local pH and the concentrations of other relevant species are illustrated in Supplementary Fig. 1. Finally, the updated local concentrations are iteratively fed back into the MKM until the variations in both pH and H2 concentration between successive iterations fall below 10−3, ensuring full convergence of the simulation. Notably, the well-established mechanisms of HER allow us to elucidate the intricate interplay between kinetics and mass transport. Although extending our established framework to electrochemical CO2 reduction reaction (eCO2RR) and electrochemical N2 reduction reaction (eNRR) is possible, it will pose additional challenges due to the complex reaction mechanisms and microenvironments at the electrodes42 as well as the complex compositions in the electrolyte.

Reaction mechanisms and free energy diagram of HER on the Au(111) surface
In this work, we utilized the Au(111) surface as a model system for our multiscale simulations of the HER mechanism due to the availability of comprehensive experimental data41 for direct comparison. We initiated our mechanistic analysis with Au(111) for two key reasons. First, the high stability and uniform surface characteristics of Au closely align with idealized DFT computational models43–45, thereby significantly enhancing the reliability of direct comparisons between theoretical simulations and experimental measurements on the Au(111) surface. Second, Au has a relatively weaker binding strength with hydrogen and exhibits no detectable underpotential hydrogen deposition (HUPD) before reaching the hydrogen reversible potential46–55, which would provide additional mechanistic insights for HER.

As shown in Fig. 2a, b (Supplementary Data 1), we plotted the free energy diagrams for the HER via both Volmer-Heyrovsky and Volmer-Tafel pathways on the Au(111) surface. The calculations were performed at 0 V vs. the standard hydrogen electrode (SHE) for the acidic pathway (H3O+ as the proton source) and −0.826 V vs. SHE for the alkaline pathway (H2O as the proton source). Note that our model distinguishes between the acidic and alkaline HER by using different explicit water structures. The proton sources are H3O+ for the acidic pathway and H2O for the alkaline pathway. Based on the significance of coverage effect in MKM, we computed the reaction energetics under both low hydrogen coverage (θ*H=0 ML) and high hydrogen coverage (θ*H=7/9 ML) conditions by varying the number of adsorbed H atoms on the surface. This comparative approach captures coverage-dependent trends across relevant electrochemical potentials. (All the computational details are summarized in Supplementary Figs. 3-8, Supplementary Data 1).

As shown in Fig. 2a, both the Volmer and Heyrovsky steps exhibit higher energy barriers in alkaline conditions than in acidic conditions. This energetic difference rationalizes the sluggish HER kinetics observed experimentally under alkaline conditions56–58. When the effect of surface coverage is taken into account, we find that higher hydrogen coverage leads to increased energy barriers for both steps, primarily due to the weakened binding strength of hydrogen. This facilitates the reverse Volmer reaction while promoting the forward Heyrovsky step. As governed by both pH conditions and surface coverage, a similar trend was also found for the Volmer-Tafel pathway in Fig. 2b. Meanwhile, the dynamic processes of the Volmer, Heyrovsky, and Tafel pathways under both alkaline and acidic conditions are shown in Fig. 2c, d. Under alkaline conditions, the explicit solvation structures are represented by a four-water cluster (H8O4) with H2O serving as the proton source, where the adjacent H2O molecules would transfer a proton to the H2O molecule that directly serves as the proton source in the Volmer and Heyrovsky steps. Under the acidic conditions, an Eigen cation (H9O4+) with a hydronium ion works as the proton source.

Multiscale simulations of polarization curves of HER on the Au(111) surface
By incorporating the complete free energy profile of the HER network into our multiscale simulation framework shown in Fig. 1b, we obtained the simulated polarization curve of HER on the Au(111) surface in Fig. 3a, which exhibits a good agreement with the experimentally measured polarization behavior in Fig. 3b41. A discrepancy is observed at pH=4 (Fig. 3b) where the simulated current slightly exceeds experimental values. While 50 mV s−1 is a standard experimental condition, we suggest that in dilute acid environments, the finite scan rate might potentially hinder complete diffusion equilibrium, possibly contributing to the observed difference. To evaluate further the robustness and reproducibility of the established multiscale simulation framework, we integrated the DFT data from Liu and coworkers39 into our mass transport model. Then, we specifically simulated the polarization curves of HER on the Pt(111) surface, which also showed good agreement with experimental results (Supplementary Fig. 9). Indeed, all the above agreements between theory and experiment suggest the robustness of our approach.

Understanding the effects of pH on polarization curves
Although our multiscale simulations can appropriately reproduce the experimental polarization curves, they only reflect the macroscopic outcome of the interplay between elementary steps and local environmental effects. Thus, the specific contribution of each factor remains unclear. Taking pH effect as an example, the HER activity in Fig. 3b exhibits strong pH dependence at acidic conditions (pH=1–4), where varying pH significantly alters the reaction kinetics. In contrast, the polarization behavior remains largely unaffected under alkaline conditions (pH=5–13).

To elucidate this pH-dependent behavior, we simulated the relationship between local pH at the cathode surface and reaction rates across the whole pH range. As shown in Fig. 3c, there are three distinct regimes, i.e., a plateau region where pH(local) ≈ pH(bulk), a pH-jump regime with proton depletion, and a linear regime exhibiting gradual pH increase with current density. Meanwhile, Fig. 3c can provide mechanistic explanations of the distinct trends observed in the polarization curves in Fig. 3a.

More specifically, at pH(bulk)=3 in Fig. 3c, the reaction proceeds slowly at low current densities (10−6 to 100 mA cm−2), allowing sufficient time for proton transfer from the bulk solution to the electrode surface. Consequently, the local environment maintains equilibrium with the bulk environment (pH(local) ≈ pH(bulk)). However, when the current density reaches 3.1 mA cm−2, the reaction rate exceeds the proton diffusion rate, resulting in progressive proton depletion near the electrode surface. This triggers a rapid transition of the pH(local) from acidic (pH=3) to alkaline (pH≈11). Notably, this pH jump is accompanied by the appearance of a distinct current plateau in the polarization curve (where current density remains constant despite potential changes). This has been widely observed experimentally in mild acidic media, which were dominated by mass transport limitations41. Then, at higher current densities in Fig. 3c, pH(local) exhibits a linear dependence on the applied current density, ultimately converging to identical trajectories regardless of pH(bulk). The coincident behavior of the polarization curves across a broad bulk pH range (5–13) can be attributed to the combined influence of two key conditions established in this regime: the convergence of the local pH to a universally alkaline regime once currents exceed 10−1 mA cm−2, and the remarkable insensitivity of the HER kinetics to its precise value within that range59. Consequently, the corresponding HER polarization curves exhibit a high degree of consistency across the investigated potential window. Our systematic analysis demonstrates that the local environment becomes dramatically different from bulk conditions (from pH=3 to pH=14) as current density increases, which will significantly impact the microkinetic modeling results if the mass transport phenomenon hasn’t been considered.

Understanding mechanisms behind the variation trend of the polarization curves under acidic conditions
The polarization curves in Fig. 3b under acidic conditions consistently exhibit a characteristic current plateau where the current density remains constant despite increasing negative potential. Notably, Nagel and Wendler60 first reported this phenomenon, while subsequent experimental and theoretical studies have consistently observed this behavior and attributed it to proton transfer limitations38–41. However, the underlying reaction mechanisms behind this phenomenon remain poorly understood. Indeed, our multiscale simulations indicate that the current plateau can be attributed to a dynamic equilibrium between the potential-dependent decrease in the energy barrier of the Volmer step and the increase in local pH at increasingly negative potentials.

As shown in Fig. 4a, we calculated the potential-dependent energy barriers (Ga) of Volmer steps in both acidic (Volmer(acid)) and alkaline (Volmer(alkaline)) conditions at high H* coverage and pH(bulk)=3. As a proton-coupled electron transfer (PCET) process, the Volmer step exhibits a remarkable decrease in the energy barrier with the potential shifting negatively, and the local pH progresses through three distinct stages (Fig. 4a). Therefore, the competition between the decrease in the energy barrier and the increase in local pH (decreased proton concentration) creates the characteristic three-stage polarization profile, featuring the observed current plateau.

More specifically, the local pH remains constant from 0 to −0.7 V vs SHE in Fig. 4a. In this respect, the continuous decrease in the energy barrier of the Volmer step with increasing potentials results in a substantial enhancement in current density. Starting from −0.7 V vs SHE, further negatively shifting the potential induces a rapid increase in local pH. However, the acidic Volmer step still serves as the sole source of H* (Fig. 4b). Although the decreasing proton concentration (reflected in an increased local pH) impedes reaction kinetics, this effect is precisely counterbalanced by the ongoing reduction in energy barrier for the acidic Volmer step, thereby resulting in a dynamic steady state that maintains a constant current density despite increasing potential. At potentials above −1.4 V vs SHE in Fig. 4b, the increase in local pH becomes gentle after the sharp jump from acid to alkaline. Then, the contribution of the acidic Volmer step decreases rapidly under such severely proton-depleted conditions. When the potential is higher than −1.8 V vs SHE, the alkaline Volmer step emerges as the exclusive mechanism for H* formation.

Furthermore, we investigated the variation of the rate-determining step across different potential ranges by examining the dependence of Tafel slopes on local pH and applied potential at bulk pH=3. As shown in Fig. 4c, the reaction kinetics are governed by three distinct factors during this process. At the low potential range, our simulated Tafel slope of 69.4 mV dec−1 is very close to the experimental value of 68 mV dec−1 as measured by Kahyarian and coworkers using Au RDE in 0.1 M NaClO4 electrolyte49, thereby confirming the Tafel step as the rate-determining step under these conditions46. In the current plateau region from −0.7 V to −1.4 V vs SHE, the Tafel slope approaches ∞ mV dec−1, indicating that the reaction is no longer limited by any PCET step but by proton transfer limitations as previously analyzed. At the high potential regime, the simulated Tafel slope of 141.5 mV dec−1 also agrees well with the experimental value of 141 mV dec−1 obtained in 0.1 M H2SO4 solution, demonstrating that the reaction becomes limited by the alkaline Volmer step52. In this respect, the water dissociation process under alkaline conditions becomes the predominant pathway for H* formation in this potential regime, as shown in Fig. 4b. Notably, between the low-potential region and the vertical plateau, a transition zone exists that exhibits intermediate Tafel slopes (ranging from 69.4 to ∞ mV dec−1), which can be attributed to mixed control by both electrochemical reaction kinetics and proton transfer. This mixed control behavior is similarly observed in the region between the vertical plateau and the high-potential Tafel line.

Building upon our systematic investigation of the potential-dependent behavior at a bulk pH of 3, which is a representative acidic condition, we have elucidated the underlying mechanistic origins of the distinct polarization curve trends. Through detailed Tafel slope analysis in Fig. 4c, we identify five characteristic potential regimes in the acidic polarization curves, each governed by different rate-limiting factors, as comprehensively mapped in Fig. 4d.

More specifically, in the low-potential regime with pH(local) ≈ pH(bulk), the reaction kinetics are exclusively limited by the sluggish Tafel step and yield negligible current density. As the potential shifts negatively, a measurable current density emerges under the mixed control of both PCET kinetics and mass transport effects, as the electrode surface environment begins to diverge from the bulk solution composition. As the potential shifts further negatively, the polarization curve exhibits a characteristic plateau, limited only by mass transport, where increasing the potential fails to enhance the current density. Upon entering the high-potential regime, the local environment stabilizes in an alkaline state distinct from the acidic bulk solution, and the local pH gradually increases with potential. As illustrated in Fig. 4b, this regime features competitive interplay between alkaline and acidic Volmer pathways, exhibiting mixed control by both PCET kinetics and mass transport, analogous to previous transitional regions. Ultimately, with a further increase in negative potential and local pH, the alkaline Volmer step emerges as the exclusive pathway for H* formation, marking the transition to a final regime governed purely by PCET kinetics. Apparently, except when the electrocatalytic process is solely controlled by the PCET process (as indicated by the red and blue dashed lines in Fig. 4d), a comprehensive understanding of catalytic processes necessitates explicit consideration of mass transport effects in theoretical modeling.

Notably, the choice of pH=3 for our detailed analysis of polarization behavior under acidic conditions was guided by fundamental electrochemical considerations across the pH range. In strongly acidic media, the high bulk concentration of H3O+ leads to rapid mass transport, resulting in polarization curves dominated by proton-coupled electron transfer kinetics without a well-defined diffusion-limited current plateau. In alkaline and neutral environments, the limited availability of H3O+ shifts the proton source to water molecules, maintaining the reaction under kinetic control throughout the potential range. In contrast, dilute acid conditions at pH=3 provide a unique window where different controlling regimes are prominently exhibited within experimentally accessible current densities. This behavior is well-documented in systematic studies of HER polarization curves61. Our multiscale framework builds upon these experimental observations by providing atomistic-scale insights into the mechanisms governing these transitions, making pH=3 an ideal case for demonstrating how our model deconvolutes the origins of different polarization curve segments.

Fine-tuning and optimizations of mass transport
Till now, we have systematically examined and illustrated the physical origins of the pH-dependent polarization behaviors of HER on the Au(111) surface. However, experimental results in Fig. 5a reveal that varying catalysts significantly alter the polarization behavior in both PCET and mixed-control regions, but do not influence the diffusion-limited current density. Therefore, we further investigated how to alleviate the diffusion limitations via tuning mass transport. As clearly shown in Fig. 5b, the increase in rotation rate will profoundly improve the diffusion-limited current density. Meanwhile, the simulated diffusion-limited current densities in Fig. 5c exhibit a linear dependence on the square root of rotation rate (ω1/2), in agreement with Levich equation predictions62.

To further elucidate how rotation influences mass transport and consequently affects reaction kinetics, we analyzed the relationship between rotation rate and both local pH and convective flux in the boundary layer (pH(bulk)=3, j=3.5 mA cm−2). As shown in Fig. 5d, our simulated results indicate that the boundary layer maintains an alkaline local pH environment at lower rotation rates (ω=400-1600 rpm), which are insufficient to adequately alleviate the diffusion limitations. However, when the rotation rate exceeds 1600 rpm, the local pH progressively decreases to 3–4, approaching the bulk pH. Notably, a further increase in the rotation rate from 2500 to 6400 rpm does not produce significant additional pH variation. Thus, our simulations indicate that 2500 rpm essentially represents the limit for rotation-enhanced mass transport. Further analysis of the convective flux (Supplementary Fig. 10) reveal that while the convective flux increases with distance from the electrode surface at a constant rotation rate, enhanced rotation does not directly affect mass transport at the electrode interface. Instead, it improves mass transport at the boundary layer, thereby indirectly modifying the interfacial reaction environment. All these findings collectively demonstrate that rotation rates within a specific range (ω=1600–2500 rpm) can optimize local reaction conditions to promote the reaction. However, since rotation only affects the convective component of mass transport, its ideal performance is ultimately limited to achieving pH(local) = pH(bulk), with no capacity for further enhancement of local pH.

Methods
Details of potential-dependent DFT calculations
All DFT calculations were performed using the Vienna ab initio simulation package (VASP)63 as implemented in the Atomic Simulation Environment (ASE)64. The Bayesian error estimation functional was selected with van der Waals correction (BEEF-vdW)65 for exchange-correlation effect. The projector augmented wave (PAW)66 method was employed to treat electron-ion interactions, with a plane-wave basis set and an energy cutoff of 400 eV for geometry optimization. A 3-layer p(3×3) supercell was employed to model the Au(111) surface, and the bottom two atomic layers of the model were fixed to mimic the bulk catalyst. The Au(111) unit cell with an area of 6.90×10−19 m2 contains 9 active sites, yielding an active site concentration of 2.16×10−5 mol m−2.

For the electrolyte, a continuum water model was applied to the bulk phase, with the system neutralized using the linear Poisson-Boltzmann implicit solvation model in VASPsol67. The electrolyte was modeled with a relative permittivity of ϵb=78.4 and a Debye length of λD=3.0 Å (benchmarked as in Supplementary Fig. 11). A 15 Å vacuum space was introduced along the z-axis to eliminate periodic image interactions. To better describe the effect of hydrogen bonds on the adsorbates, our model combines explicit and implicit water molecules. The hydronium ion is explicitly solvated by three water molecules to form an eigen cation (H9O4+), thereby modeling acidic conditions through H3O+ as the proton source, while alkaline conditions are modeled using a four-water cluster (H8O4) with H2O as the proton source. The hybrid explicit/implicit solvation model builds upon established computational electrochemistry approaches67–69, combining implicit electrolyte screening with explicit water molecules to capture essential hydrogen bonding effects. The methodology is supported by previous benchmarks confirming its reliability for interfacial simulations70,71, providing an ideal balance between accuracy and efficiency for investigating multiscale electrochemical phenomena. While these static calculations cannot fully capture the dynamic hydrogen bond fluctuations and presolvation effects highlighted in recent studies72, the modest energetic variations suggest that our primary conclusions regarding multiscale framework validation and key polarization behaviors analysis remain valid. We emphasize that incorporating such dynamic solvation effects represents an important direction for future methodological development in electrochemical interface modeling. Structural optimizations converged at forces below 0.05 eV Å−1 and energy changes below 10−4 eV, with Brillouin zone sampling using a 2×2×1 Monkhorst-Pack k-point mesh. Transition states were located using the climbing-image nudged elastic band (CI-NEB) method73 with 8 images. Vibrational frequency analysis was employed to calculate zero-point energies (ZPE) and entropic contributions to the free energy at 298.15 K, considering exclusively the 3N degrees of freedom associated with the directly reacting adsorbates.

For the constant electrode potential (CEP) modeling, we used the fully converged constant-potential (FCP) algorithm to perform grand canonical ensemble modeling with density functional theory (DFT) calculations as developed by Xiao74. The system’s electron count was modulated to control the Fermi energy, maintaining both a fixed work function and electrode potential to achieve the target state. The reductive potential U is related to the work function (Wf) of the system and absolute standard hydrogen electrode potential (ϕSHE) by,

U(vs. SHE) = Wf / e − ϕSHE (1)

the standard hydrogen electrode potential (ϕSHE) was set to 4.43 V75. In electrochemical systems, the free energy is conventionally expressed as a function of the applied electrode potential,

G(U) = EDFT + ZPE − TS − qμe (2)

the grand canonical character of the electrochemical system is embedded in the final term through the net charge per supercell (q), reflecting electron density fluctuations and the μe is the chemical potential under the target potential. Here, EDFT, ZPE, and TS correspond to the DFT electronic energy, vibrational zero-point energy, and entropic effects, respectively. To bridge constant-charge and constant-potential free energies, we performed quadratic potential-dependent fitting for all reaction species, which will be shown in Supplementary Tables 5 and 6. This approach quantitatively relates the applied potential to both the reaction free energy (ΔG(U)) and the activation energy barrier (Ga(U)).

Free energy diagrams
The free energy diagrams (Fig. 2a, b) use the initial Volmer step state as energy zero (0 eV): For the acidic pathway, the reference is * + H3O+ + e−, and the alkaline pathway, the reference is * + H2O + e−.
All subsequent states are calculated relative to this reference. To evaluate coverage effects, we compared two surface environments: θ*H=0 ML for the clean surface without pre-adsorbed H, and the θ*H=7/9 ML surface with seven pre-adsorbed H atoms in a unit cell. Corresponding structures are shown in Supplementary Figs. 5-6.

Microkinetic modeling
Mean-field microkinetic modeling was performed by incorporating potential-dependent DFT data, local pH, and H2 concentration at the cathode surface to determine reaction rates for all elementary steps and surface coverages of reaction species. The local pH is defined as the pH at the cathode surface, which is directly calculated from the local concentration of hydronium ions (cH3O+,local) at that plane, schematically illustrated in Supplementary Fig. 1. The microkinetic modeling utilized energetics calculated at 0 ML hydrogen coverage as reference, with linear interpolation between 0 ML and 7/9 ML datasets to account for coverage effects. The model considers reactant and product adsorption/desorption processes along with all elementary reaction steps, all possible reaction mechanisms are depicted below,

H3O+ + * + e− ⇌k ±1 H* + H2O (3)

H3O+ + H* + e− ⇌k ±2 H2(aq) + H2O + * (4)

H2O + * + e− ⇌k ±3 H* + OH− (5)

H2O + H* + e− ⇌k ±4 H2(aq) + OH− + * (6)

H* + H* ⇌k ±5 H2(aq) + 2* (7)

Transition state theory (TST) was employed to compute reaction rates (ki) for each elementary step i,

ki = kB T / h * e −Ga,i(U) / kBT (8)

where kB is the Boltzmann constant, T is the temperature of the system, h is the Planck constant, and Ga,i(U) is the activation energy barrier at the applied potential. The equilibrium constant Ki of each elementary reaction i will be described by,

Ki = e −ΔGr,i(U) / kBT (9)

where ΔGr,i(U) is the reaction free energy of the elementary reaction i, which is evaluated from the GC-DFT calculation. The model of the reaction network is usually related to the steady-state of a set of coupled nonlinear ordinary differential equations (ODE),

∂θi / ∂t = 0, Σθi = 1 (10)

where θi is the coverage of species i on the catalyst, The standard concentration cΘi = 1 mol L−1 was applied to all species i, while water was treated at its bulk reference concentration cΘH2O = 55.6 mol L−1. The coverage will be described by the following equations:

∂θH / ∂t = k1 θ* cH3O+ − k−1 θH cH2O + k2 θH cH3O+ − k−2 θ* cH2O cH2 + k3 θ* cH2O − k−3 θH cOH− − k4 θH cH2O + k−4 θ* cOH− cH2 − 2k5 θH2 + 2k−5 θ*2 cH2 (23)

and,

θ* + θH = 1 (24)

Definition and calculation of volmer step ratio
The Volmer step ratio quantifies the relative contribution of the acidic versus alkaline proton transfer pathways to the formation of the H* intermediate. It is defined as the percentage of the total H* production rate that proceeds through the acidic Volmer mechanism:

Volmer step ratio = rVolmer(acid) / (rVolmer(acid) + rVolmer(alkaline)) × 100% (25)

The individual reaction rates are obtained from the solution of the microkinetic model, which is self-consistently coupled with the continuum transport calculations. The potential dependence of this ratio reflects the shifting availability of H3O+ versus H2O as proton sources under different interfacial environments.

Continuum transport model
A schematic illustration of the one-dimensional rotating disk electrode (RDE) system is presented in Supplementary Fig. 1. The Au(111) cathode is used for the reduction of H2O to H2, and the thickness of the boundary layer (δ) is equal to 4×10−4 m. Species conservation is invoked to solve for concentrations for five of the species (H2(aq), H2O, H3O+, K+, ClO4−) in the electrolyte within the mass transport domain, all the species satisfy the conservation of:

∇·Ni = Ri (26)

where Ni and Ri are the flux and the net reaction rate in the solution of species i, respectively. The total molar species flux is governed by the Nernst-Planck equation under dilute solution theory, incorporating contributions from diffusion, migration, and convection76,

Ni = −Di ∇ci − zi ui ci ∇ϕ + vci (27)

the three terms on the right side of the equation quantitatively describe the diffusion, migration, and convection phenomena, respectively. Where Di, ci, and zi represent the diffusion coefficient, concentration, and charge of species i, respectively, while ϕ and x denote the electrolyte ionic potential and position. The mobility can be determined from the Nernst-Einstein relationship:

ui = Di / RT (28)

The migration term applies to charged ions: H3O+, OH−, K+, ClO4−. F is the Faraday constant, R is the ideal gas constant, and T is the temperature. The fluid velocity v in the RDE system is described by:

v = −0.51 ω3/2 υ−1/2 x2 (29)

where ω is the rotation speed and υ is the kinematic viscosity of the solution, which are maintained at 1600 rpm and the experimental conditions. The electrolyte current density jl can be obtained from the ionic flux of all ions:

jl = F Σ zi Ni (30)

and all ionic profiles must satisfy the electroneutrality condition,

Σ zi ci = 0 (31)

Furthermore, the bulk ionization of H2O should be considered in our model:

2H2O ⇌k ±w H3O+ + OH− (32)

where k ±w represent the forward/reverse reaction rates of water dissociation, the value of the equilibrium constant Kw = 1.0 × 10−14, the forward reaction rate constant for water ionization k+w = 2.4 × 10−2 mol m−3 s−1 37.

The cell potential applied to the RDE comprises the thermodynamic potential, kinetic overpotential, ohmic overpotential, and Nernstian overpotential. The ohmic losses can be described by:

ηohmic = ∫ δ0 (j / κ) dx (33)

where κ is electrolyte conductivity. The deviation of the solution concentration at the catalyst surface from the bulk solution causes the Nernstian overpotential and results in an increase in the equilibrium potentials. The Nernstian overpotential at the cathode is:

ηner,HER = 2.303 RT / F (pHlocal − pHbulk) − RT / 2F ln(cH2,local / cH2,bulk) (34)

where cH2,local and cH2,bulk are the concentrations of H2 in the bulk electrolyte and at the cathode, respectively. The Nernstian overpotential of HER at the cathode comprises pH loss and H2 production, which result from the rapid change of pH and H2 concentration at the cathode.

Multiscale simulations
Supplementary Fig. 2 presents our multiscale simulation framework, which integrates the continuum transport model, microkinetic model, and GC-DFT. The simulation protocol initiates with the applied cathode potential (vs. SHE) as the primary input. GC-DFT calculations yield key kinetic and thermodynamic parameters, including activation barriers and reaction free energies for elementary steps. These parameters directly inform the rate constants (k), equilibrium constants (K), coverage (θ), and reaction rate (r) for all considered surface processes. The concentration of H2 (cH2) and pH are initialized using their bulk values. The MKM parameterized with GC-DFT-derived energetics, then computes the potential-dependent partial current densities for all products. These current densities are finally fed into the continuum transport model, which solves the coupled mass transport effects to obtain the local concentrations of H2 and pH at the cathode-electrolyte interface. This MKM-transport coupling typically converges in fewer than four iterations in strongly acidic/alkaline conditions, but requires at least 20 iterations in dilute acid conditions (pH 3-5) due to the sharp local pH changes near the diffusion-limited region. Convergence was achieved when the relative change in local pH and H2 concentration fell below 10−3.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors developed a multiscale simulation framework integrating grand canonical density functional theory (GC-DFT), microkinetic modeling, and continuum transport analysis to reproduce experimental hydrogen evolution reaction polarization curves across all pH ranges on Pt(111) and Au(111), thereby revealing potential-dependent local environment evolution and explaining how mass transport influences the formation of characteristic current plateaus.\n\nComment 1. How do the authors define “local pH”? Could the main text be supplemented with a precise mathematical formula for this definition? Adding such details would enhance clarity and broaden the article's accessibility to readers.\nComment 2. The sentence “all the details of our framework……to applications in electrochemical CO2 and N2 reduction.” Appeared nearly identically twice. While the manuscript highlights the method's generalizability, potential challenges may arise in its application to electrochemical CO2 and N2 reduction, given the more complex local environments in these systems.\nComment 3. On page 17, section “Details of Potential-dependent DFT Calculations.” the authors use a relative permittivity of ϵb = 78.4 and a Debye length of λD = 3.0 Å. However, under constant potential, the solution’s relative permittivity near the electrode surface is typically much lower than 78.4 [see J. Phys. Chem. Lett. 14, 685−693 (2023) and J. Chem. Phys. 19, 391–394 (1951)], and local pH variations may influence the Debye length. Could these fixed parameters introduce significant errors in the calculation results?\nComment 4. In Figure 3B, the current plateau at pH = 4 is noticeably higher than in Figure 3A. What could be the underlying cause of this discrepancy?\nComment 5. They employ Eigen cation with a hydronium ion as the proton source. However, the results of thermodynamically-based DFT calculations do not adequately account for kinetic processes at room temperature. Are factors such as fluctuations of hydrogen bond networks and system charge significant? For example, as noted in the reference J. Phys. Chem. Lett. 2025, 16, 4572−4578 (“Presolvation Dynamics Preceding the Hydrated Proton Transfer in the Electrical Double Layer”), while the Eigen structure is the primary form of protonated species under acidic conditions, a partial transition from Eigen (H2O)3H3O+ to pentamer (H2O)4H3O+ occurs during proton transfer.\nComment 6. The reaction energy of Volmer-Tafel pathway under low hydrogen coverage (θ∗H=0 ML) and high hydrogen coverage (θ∗H=7/9 ML) should be equal in Fig 2b.\nComment 7. The free energy of proton is very sensitive to the number of explicit waters used in the model. The water number used in the work is not enough to obtained stable proton, which bring errors in the thermodynamic energy calculation. One consequence is that the reaction energy of the whole pathway is negative even the author claimed the reaction is operated under equilibrium potential.", "reviewer_2": "The work by Shou et al. simulates the hydrogen evolution reaction primarily on a Au electrode using a multiscale approach that combines grand canonical density functional theory (GC-DFT), the microkinetic model (MKM), and a continuum transport model. The framework looks promising and seems well-suited for a wide range of applications, as the authors suggest. However, the manuscript requires significant revision before publication.\n1. It's unclear how many iterations between the MKM and the transport model are necessary to achieve converged results, based on the framework the authors have proposed. The authors calculated two different hydrogen coverages (0 and 7/9 ML), but it's not clear which results were used for the subsequent MKM simulation. It appears the authors may have ignored the effect of changing hydrogen coverage, as there is no explanation for how these specific coverages were chosen. At a minimum, the authors should calculate the surface Pourbaix diagram to acknowledge this important interface phenomenon, even if it doesn't fully capture every detail.\n2. The hybrid explicit/implicit model used for the GC-DFT calculation lacks supporting evidence. The authors need to discuss the origin of this model and provide a rigorous test to demonstrate its validity. This is crucial for the scientific community to trust the model as an accurate methodology for simulating solid-liquid interfaces. The current approach appears to be heavily biased by the authors' perspective.\n3. The statement on page 10, \" This convergence explains why polarization curves exhibit nearly identical trends at pH(local) = 5–13 when currents exceed 10−1 mA/cm2, as…\", is confusing. Based on Figure 3c, I observe a deviation from this explanation. At this current density, the local pH values for bulk pH 5 and 13 appear different (e.g., local pH 9 for bulk pH 5 and local pH 13 for bulk pH 13).\n4. In the section titled \" Understanding the Mechanism of Polarization Curve Variation Trend under Acidic Conditions,\" the authors discuss results by setting the bulk pH to 3. It would be more impactful to show results for a strong acid, a strong base, and a neutral solution (e.g., pH 0, 14, and 7). The authors also need to rationalize why they chose this specific pH value.\n5. Although the authors discuss their methodology, they do not provide any input/output files needed to reproduce their results. To ensure reproducibility and transparency, they must share the raw input and results data for the GC-DFT, MKM, and transport modeling on a public platform like a GitHub page. This is essential and should be shared during the review process.", "reviewer_3": "The paper “Toward Rational Understanding of the Hydrogen Evolution Polarization Curves Through Integrating Ab Initio Simulations and Mass Transport Effects” presents a valuable modeling study aimed at understanding polarization curves for the hydrogen evolution reaction (HER) on Au(111) and Pt(111) surfaces. The authors employ multiscale simulations that integrate grand-canonical density functional theory (GC-DFT) calculations, a microkinetic model, and a continuum mass transport model.\n\nThe central contribution of the paper is the demonstration that mass transport phenomena significantly influence polarization behavior and catalyst performance, and that their explicit inclusion in simulations is essential to accurately reproduce experimental polarization curves.\n\nOverall, the work addresses a timely and important topic in the theoretical description of electrochemical processes. The manuscript contains sufficient novelty to merit consideration for publication in Nature Communications. However, several issues should be addressed prior to acceptance:\nIn the Abstract and Introduction, the authors state that results are presented for both Au(111) and Pt(111). However, the Results section primarily discusses Au(111). I recommend either including Pt(111) data in the main text or removing its mention from the Abstract.\nMany sentences in the Abstract and Introduction are repetitive. Please revise to improve conciseness and avoid redundancy.\nThe reference energy of the products used in the free energy diagrams (Fig. 2a and Fig. 2b) is not clearly defined. The method of determining these reference energies should be explicitly stated.\nThe definition and calculation procedure of the Volmer step ratio (%) reported in Fig. 4b are unclear and require clarification.\nThe manuscript introduces the concept of “locality” of the pH. How is this defined? If it corresponds to the concentration of H₃O⁺ within the boundary layer, please provide the boundary layer thickness, at least for the case of an RDE rotation speed of 1600 rpm.\nThe flow of argumentation in the Results section is at times difficult to follow. Improving the clarity and organization of this section would strengthen the paper.\nThe expression for the time derivative of H* coverage (Eq. 15) is formally incorrect. The two terms associated with the Tafel step should be multiplied by a factor of 2 (see Eq. 115 in [Int. J. Hydrogen Energy, 2019, DOI: 10.1016/j.ijhydene.2019.05.183]). Please verify whether this is a typographical error or requires correction.\nThe manuscript reports Eq. 21 as j(H₂) = 2·F·r(H₂), with r(H₂) = r₂ + r₄ + r₅. However, r(H₂) is originally defined as r(H₂) = r₁ + r₂ + r₃ + r₄, which can be reformulated as the authors present only under the steady-state assumption (dθ(H*)/dt = r₁ – r₂ + r₃ – r₄ – 2·r₅ = 0). The authors should clarify whether the steady-state approximation has been explicitly applied in the microkinetic model."}, "author_response": "General Remarks: We are deeply grateful to all the reviewers for their professional and constructive suggestions, which dramatically improved the quality of our work. We have carefully addressed all the comments and revised the manuscript accordingly. The revisions are highlighted in the manuscript. The point-by-point responses are listed below, with the reviewers’ original comments in black, our responses in blue, and corresponding revisions of the manuscript in yellow.\nReviewer #1 (Remarks to the Author):\nThe authors developed a multiscale simulation framework integrating grand canonical density functional theory (GC-DFT), microkinetic modeling, and continuum transport analysis to reproduce experimental hydrogen evolution reaction polarization curves across all pH ranges on Pt(111) and Au(111), thereby revealing potential-dependent local environment evolution and explaining how mass transport influences the formation of characteristic current plateaus.\nAuthor Reply: We thank the reviewer for the professional comments and suggestions on our work, which have significantly improved its quality. Accordingly, we have revised our manuscript in response to your valuable suggestions. We sincerely hope our revision will address all your valuable concerns.\nComment 1. How do the authors define “local pH”? Could the main text be supplemented with a precise mathematical formula for this definition? Adding such details would enhance clarity and broaden the article's accessibility to readers.\nAuthor Reply: We thank the reviewer for this great suggestion. Indeed, the definition of the \"local pH\" is critical in our work. In this study, the local pH is defined as the pH value at the cathode surface. It is directly calculated from the local concentration of hydronium ions (c ) at that plane using the standard formula:\npH = −log (c )\n(local) 10 H3O+,local\nWhere the unit of c is mol/L. As schematically illustrated in Supplementary Figure. 1, c denotes the concentration of H₃O⁺ ions at the cathode surface. As an output of our continuum transport model, the local concentration of H₃O⁺ is determined by solving the Nernst-Planck equation (detailed in the Methods section) under steady-state conditions. Then the value of c at the position of the cathode surface is taken as c , which is then converted to pH . In essence, the pH is an iterative result obtained by simultaneously considering the dynamic equilibrium between surface reaction kinetics and mass transport of species in the electrolyte. It is therefore dependent on the applied potential, the bulk pH, the reaction rate, and the hydrodynamic conditions. We have added the specific definition of local pH in our revised manuscript.\nPlease see the lines 83-86 on page S4 in the revised SI and also below:\n“The local pH is defined as the pH at the cathode surface, which is directly calculated from the local concentration of hydronium ions (c ) at that plane schematically illustrated in Supplementary Figure 1.”\nComment 2. The sentence “all the details of our framework……to applications in electrochemical CO and N reduction.” Appeared nearly identically twice. While the manuscript highlights the method's generalizability, potential challenges may arise in its application to electrochemical CO and N reduction, given the more complex local environments in these systems.\nAuthor Reply: We sincerely thank the reviewer for bringing this crucial point to our attention. We also apologize for the unintentional duplication of the sentence, which has been removed from the revised manuscript. Meanwhile, we fully agree with the reviewer that extending our framework to more complex reactions, such as CO and N reduction, will pose additional challenges. We have toned down our claim and discussed potential challenges in the revised manuscript.\nPlease see the lines 124-128 on page 6 in the revised manuscript and also below:\n“Notably, the well-established mechanisms of HER allow us to elucidate the intricate interplay between kinetics and mass transport. Although extending our established framework to eCO RR and eNRR is possible, it will pose additional challenges due to the complex reaction mechanisms and microenvironments at the electrodes42 as well as the complex compositions in the electrolyte.”\nComment 3. On page 17, section “Details of Potential-dependent DFT Calculations.” the authors use a relative permittivity of ϵb = 78.4 and a Debye length of λ = 3.0 Å. However, under constant potential, the solution’s relative permittivity near the electrode surface is typically much lower than 78.4 [see J. Phys. Chem. Lett. 14, 685−693 (2023) and J. Chem. Phys. 19, 391–394 (1951)], and local pH variations may influence the Debye length. Could these fixed parameters introduce significant errors in the calculation results?\nAuthor Reply: We thank the reviewer for this professional comment. We also agree that the relative permittivity near the electrode surface is typically lower than that in bulk. Meanwhile, we investigated the effects of varying ϵ on the calculated results (supplementary Figure 11). Although the conclusion is that there are no significant errors in the calculated results, we would say that such comparisons are inadequate due to the high complexity of the interfacial dielectric environment. Therefore, we hope to share with you our understanding and approach to tackling this challenging issue.\nIn our work, a continuum water model was applied to the bulk phase, and the system was neutralized using the linear Poisson–Boltzmann implicit solvation model in VASPsol. The parameters adopted in our calculations (ϵ =78.4 and λ =3.0Å) correspond to a bulk aqueous electrolyte at room temperature with approximately 1 M ionic strength. In the interface region where the electronic charge density increases rapidly from zero toward the values inside the solute, the properties of the implicit solvation model, such as the relative permittivity (ϵ ) and the Debye screening length (λ ), gradually change from their bulk electrolyte values to those of vacuum in the explicit solute region. This scaling is governed by the following shape function:\n𝜁[𝑛(𝑟⃗)] = 1/2 erfc{...}\nwhich ensures a smooth transition of electrolyte properties across the solute–electrolyte boundary. This treatment captures the interfacial decay of the dielectric environment. Directly modifying ϵ in the input would artificially change the bulk electrolyte properties. Indeed, the effective permittivity at the electrode interface is often significantly lower than the bulk value, and the current model, which relies primarily on the spatial shape function 𝜁.\nNotably, both the relative permittivity (ϵ ) and the Debye screening length (λ ) depend on local concentration of spices. Specifically, the permittivity can be expressed as, ϵ = ϵ , and the Debye length as, λ = ...\nVariations in the local ionic concentration (c) thus directly modulate both ϵ and 𝜆 , and thereby alter the screening strength near the interface. While these relationships highlight the complex coupling between ion distributions and dielectric properties, they also indicate that neglecting such concentration-dependent effects may introduce quantitative inaccuracies in modeling the interfacial environment.\nNevertheless, within our multiscale simulation framework, potential-dependent reaction energetics are first obtained from the GC-DFT and then passed into the MKM. The MKM is coupled with continuum transport equations and solved iteratively to yield local reaction environments and reaction rates. Compared with DFT calculations, solving the MKM ODEs and finite-element transport equations is computationally much more tractable, making such iterative schemes practical while still capturing the influence of the local environment on kinetics and coverages.\nTo achieve fully quantitative accuracy, one would ideally incorporate the feedback of the local environment into the GC-DFT step itself—namely, using the electrolyte concentrations obtained from the transport model to update the dielectric response and Debye length in the DFT electrolyte model, and then iterating until convergence. This strategy would allow the dielectric constant and screening length to evolve self-consistently with the local interfacial environment. However, such a fully integrated workflow is extremely challenging to realize in practice and, to our knowledge, has not yet been achieved in existing multiscale electrochemical simulations [Proc Natl Acad Sci, 2017, 114, E8812-E8821; ACS Catalysis, 2019, 9, 6194-6201; J Am Chem Soc, 2020, 142, 4985-4989; Chem Catalysis, 2021, 1, 1331-1345].\nWhile our multiscale framework does not dynamically update the relative permittivity (ϵ = 78.4) and Debye length (λ = 3.0 Å) based on the local environment, we have conducted sensitivity analyses to evaluate their impact. As highlighted in theoretical studies of CO₂RR [Energy & Environmental Science, 2019, 12, 3380-3389] , where the local relative permittivity was shown to drop to 50 under high electrode potentials. Following this insight, we systematically tested the relative permittivity of 78.4, 70, 60, and 50 with their corresponding Debye lengths (3.000, 3.175, 3.429, and 3.757 Å, respectively) for the (H O ) structure on Au(111). As shown in Figure R1, the resulting variation in reaction energetics is minor—for instance, the energy difference between ϵ = 78.4 and ϵ = 50 is 0.075 eV, which is within the typical uncertainty of DFT-based calculations.\nFigure R1 | Sensitivity analysis of dielectric constant effects on energetics. Energy profiles for the (H₈O₄) structure on Au(111) calculated with different relative permittivity (ϵ = 78.4, 70, 60, 50) and their corresponding Debye lengths (λ = 3.000, 3.175, 3.429, 3.757 Å).\nTherefore, while we acknowledge that a method capable of dynamically updating the relative permittivity and Debye length with the local environment would be ideal, the observed insensitivity of key energetics to these parameters confirms that the fixed values used in our study do not introduce significant errors or alter the main conclusions.\nWe envision that the future development of a unified, open multiscale simulation platform integrating GC-DFT, MKM, and continuum transport models will enable the precise consideration of local environmental effects on both the relative permittivity and the Debye length at electrified interfaces.\nWe have added Figure R1 as Supplementary Figure 11 in the revised SI and revised the manuscript (line 412 on Page 17) accordingly.\nComment 4. In Figure 3B, the current plateau at pH = 4 is noticeably higher than in Figure 3A. What could be the underlying cause of this discrepancy?\nAuthor Reply: We thank the reviewer for this insightful question. We acknowledge that there is a discrepancy between the simulated and experimental current plateaus at pH = 4. This discrepancy could primarily be attributed to kinetic limitations in the experimental system, which prevent it from reaching a complete steady state —a fundamental assumption in our continuum transport model.\nMore specifically, a plausible explanation is related to the scan rate used in the experiment (50 mV/s). Our model calculates the polarization curve under steady-state conditions, where the diffusion and reaction rates are fully equilibrated at each potential. However, at pH = 4, the concentration of H₃O⁺ in the bulk is very low (~10⁻⁴ M). Under a finite scan rate, the diffusion of H₃O⁺ from the bulk to the electrode surface may not be sufficiently fast to keep pace with the rapidly increasing reaction rate driven by the changing potential. This leads to a sharp rise in the consumption rate of H₃O⁺ at the surface, which can temporarily exceed the maximum supply rate achievable by diffusion from the bulk solution. Consequently, the H₃O⁺ concentration at the electrode surface in the experiment can be depleted to a level even lower than that in the steady-state model prediction. A higher local pH (lower H₃O⁺ concentration) directly leads to a lower observed diffusion-limited current, which is exactly the discrepancy seen in Figure 3B (experiment) compared to 3A (simulation). While other factors, such as the electroneutrality assumption in the transport model and the inherent limitations of the DFT description of the interface, may also contribute to minor deviations, we believe the non-steady-state effect due to the combination of low concentration and finite scan rate is the dominant cause of this specific discrepancy.\nWe have revised our manuscript to include a discussion of this limitation. We also explicitly stated that the simulation represented an ideal steady state, and that deviations at low concentrations could be expected under typical experimental scanning rates due to transient mass transport effects.\nPlease see the lines 180-184 on pages 8 in the revised manuscript and also below:\n“A discrepancy is observed at pH = 4 (Fig. 3B) where the simulated current slightly exceeds experimental values. While 50 mV/s is a standard experimental condition, we suggest that in dilute acid environments, the finite scan rate might potentially hinder complete diffusion equilibrium, possibly contributing to the observed difference.”\nComment 5. They employ Eigen cation with a hydronium ion as the proton source. However, the results of thermodynamically-based DFT calculations do not adequately account for kinetic processes at room temperature. Are factors such as fluctuations of hydrogen bond networks and system charge significant? For example, as noted in the reference J. Phys. Chem. Lett. 2025, 16, 4572−4578 (“Presolvation Dynamics Preceding the Hydrated Proton Transfer in the Electrical Double Layer”), while the Eigen structure is the primary form of protonated species under acidic conditions, a partial transition from Eigen (H2O)3H3O+ to pentamer (H2O)4H3O+ occurs during proton transfer.\nAuthor Reply: We thank the reviewer for this insightful comment and for bringing to our attention this great JPCL work by Prof. Gao and Prof. Wang, which highlights the important role of presolvation dynamics and hydrogen bond network fluctuations in proton transfer. We also appreciate their contributions to the development of fixed potential methods in modeling interfacial electrochemistry.\nTo directly address the reviewer's concern about the influence of the Eigen-to-pentamer transition on our results, we have conducted additional calculations in the spirit of the referenced study (Figure R2). More specifically, we compared the energy barriers for the Volmer step using both the Eigen cation (H O ⁺) and the pentamer cation (H O ⁺) structures as the proton source across a range of potentials. The results (Figure R2) show that the maximum difference in the activation energy barrier is ~0.05 eV. This suggests that choosing either the Eigen or pentamer structure as the initial proton donor may not alter our primary conclusions, as our study focuses on establishing a multi-scale framework and identifying qualitative trends in polarization curves. However, we fully acknowledge that these dynamic fluctuations and solvation effects are crucial for determining the precise kinetic rates and would be a critical consideration for future high-accuracy dynamic simulations. We have clarified this crucial point in our revised manuscript.\nPlease see the lines 425-431 on page 18 in the revised manuscript and also below:\n“While these static calculations cannot fully capture the dynamic hydrogen bond fluctuations and presolvation effects highlighted in recent studies68, the modest energetic variations suggest that our primary conclusions regarding multiscale framework validation and key polarization behaviors analysis remain valid. We emphasize that incorporating such dynamic solvation effects represents an important direction for future methodological development in electrochemical interface modeling.”\nFigure R2. | Comparative DFT analysis of the Volmer step using Eigen (H O ⁺) and Pentamer (H O ⁺) cation structures as the proton source. (A) Comparison of the calculated activation energy barriers (G ) for the Volmer step as a function of applied potential. (B) Atomic structure models illustrating the initial state (IS), transition state (TS), and final state (FS) for the Volmer reaction pathway using the Pentamer (H₉O₄⁺) structure. Top and side views are provided for each state.\nComment 6. The reaction energy of Volmer-Tafel pathway under low hydrogen coverage (𝜃∗𝐻=0 ML) and high hydrogen coverage (𝜃∗𝐻=7/9 ML) should be equal in Fig 2b.\nAuthor Reply: We sincerely thank the reviewer for this insightful comment, which has allowed us to identify and correct an inconsistency in the representation of the reaction energies in the original Fig. 2B.\nThe reviewer is correct in principle that the overall reaction energy for the Volmer-Tafel pathway (2(H⁺ + e⁻) → H₂) should be thermodynamically invariant with coverage. The discrepancy observed in the original figure arose from an inconsistent definition of the reference state for the Tafel step between the two coverage cases. We apologize for this oversight in our initial presentation and are grateful to the reviewer for prompting this critical clarification. Fig. 2 is updated in the revised manuscript and is also shown below.\nFig. 2 | Potential free energy profiles and reaction mechanism for HER. Grand canonical (GC) potential free energy diagrams for Volmer-Heyrovsky pathway (A) and Volmer-Tafel pathway (B) on the Au(111) surface under both clean (𝜃 =0 ML) and high-coverage (𝜃 =7/9 ML) conditions at 0 V vs. SHE (H O⁺ as proton donor) and -0.826 V vs. SHE (H O as proton donor), and the energy barriers of each step are indicated. Dynamic processes of Volmer-Heyrovsky (C) and Volmer-Tafel (D) mechanisms. The yellow, red, and white balls represent Au, oxygen, and hydrogen, respectively. Blue balls represent hydrogen involved in the reaction process.\nComment 7. The free energy of proton is very sensitive to the number of explicit waters used in the model. The water number used in the work is not enough to obtained stable proton, which bring errors in the thermodynamic energy calculation. One consequence is that the reaction energy of the whole pathway is negative even the author claimed the reaction is operated under equilibrium potential.\nAuthor Reply: We thank the reviewer for raising this critical point regarding the solvation of protons and its impact on thermodynamic calculations. We agree that the solvation free energy of a proton is highly sensitive to the explicit solvation environment, and we appreciate the reviewer's careful observation. We apologize for the unclear description of our models and approach, which may have resulted in potential misinterpretation. We hope to clarify them as follows:\n1. Primary Goal of the Free Energy Diagrams: We would like to clarify that the primary purpose of Figs. 2A and 2B is not to report the absolute reaction energy with quantitative accuracy at a specific pH, but rather to compare the thermodynamic profiles between the acidic and alkaline proton sources within a consistent computational framework. The key insight we aim to convey is the relative energy difference when the proton donor changes from H O⁺ to H O, which significantly influences the reaction landscape.\n2. Clarification on pH Conditions and Model: The reviewer is correct to point out a lack of clarity in our original manuscript. The models we used are structural representations of the local proton-transfer environment rather than exact simulations of bulk pH=0 or pH=14. The statement \"at their thermodynamic equilibrium potentials, i.e., 0 V vs. SHE for acidic pathway (pH=0) and -0.826 V vs. SHE for alkaline pathway (pH=14)\" was an oversimplification that could mislead readers. The equilibrium potential is inherently tied to the model's definition of the proton chemical potential.\n3. Addressing the “Negative Reaction Energy”: The negative overall reaction energy observed in the diagrams is a known consequence of using a limited explicit water cluster to solvate the proton, which does not fully capture its absolute solvation energy. However, since the same solvation environment is used for both the initial and final states in our comparative analysis, the systematic error is partially canceled out when examining the relative trends and barriers between the acidic and alkaline mechanisms. This approach is suitable for the qualitative comparison we intended.\nWe will clarify that the free energy diagrams are intended for a qualitative comparison of reaction pathways under different proton-source conditions, and that the absolute reaction energies are subject to uncertainty due to the finite explicit solvent model. We believe these revisions will prevent misunderstanding and better frame the objective of our DFT calculations.\nPlease see the lines 140-144 on page 7, along with the corresponding figure captions for Fig. 2 in the revised manuscript and also below:\n“As shown in Figs. 2A-2B, we plotted the free energy diagrams for the HER via both Volmer-Heyrovsky and Volmer-Tafel pathways on the Au(111) surface. The calculations were performed at 0 V vs. standard hydrogen electrode (SHE) for acidic pathway (H O⁺ as proton donor) and -0.826 V vs. SHE for alkaline pathway (H O as proton donor).”\nFig. 2 | Potential free energy profiles and reaction mechanism for HER. Grand canonical (GC) potential free energy diagrams for Volmer-Heyrovsky pathway (A) and Volmer-Tafel pathway (B) on the Au(111) surface under both clean (𝜃 =0 ML) and high-coverage (𝜃 =7/9 ML) conditions at 0 V vs. SHE (H O⁺ as proton donor) and -0.826 V vs. SHE (H O as proton donor), and the energy barriers of each step are indicated. Dynamic processes of Volmer-Heyrovsky (C) and Volmer-Tafel (D) mechanisms. The yellow, red, and white balls represent Au, oxygen, and hydrogen, respectively. Blue balls represent hydrogen involved in the reaction process.\nReviewer #2 (Remarks to the Author):\nThe work by Shou et al. simulates the hydrogen evolution reaction primarily on a Au electrode using a multiscale approach that combines grand canonical density functional theory (GC-DFT), the microkinetic model (MKM), and a continuum transport model. The framework looks promising and seems well-suited for a wide range of applications, as the authors suggest. However, the manuscript requires significant revision before publication.\nAuthor Reply: We are grateful to the reviewer for finding our work promising and providing professional and constructive suggestions to improve it. Accordingly, we have revised our manuscript in response to your valuable suggestions, which made it more impactful and clearer. We sincerely hope our revision will address all your valuable concerns.\nComment 1. It's unclear how many iterations between the MKM and the transport model are necessary to achieve converged results, based on the framework the authors have proposed. The authors calculated two different hydrogen coverages (0 and 7/9 ML), but it's not clear which results were used for the subsequent MKM simulation. It appears the authors may have ignored the effect of changing hydrogen coverage, as there is no explanation for how these specific coverages were chosen. At a minimum, the authors should calculate the surface Pourbaix diagram to acknowledge this important interface phenomenon, even if it doesn't fully capture every detail.\nAuthor Reply: We thank the reviewer for these critical and insightful comments regarding the technical implementation of our model. We apologize for the unclear descriptions of simulation details, which have caused confusion. Following your valuable suggestions, we have clarified our approach and addressed each point in our revised manuscript, and also below:\n1. Convergence of the MKM-Transport Iteration. The number of iterations required for convergence is not fixed but depends on the physicochemical conditions, specifically the coupling strength between reaction kinetics and mass transport.\nMore specifically, in strongly acidic (e.g., pH=1) or strongly alkaline (e.g., pH=14) conditions, the concentration of the primary reacting species (H₃O⁺ or H₂O) in the bulk is high. Consequently, mass transport has a negligible impact on the local environment across the simulated potential window. In these cases, the system typically converges in less than four iterations because the local environment remains largely unchanged from the initial guess.\nHowever, in dilute acid conditions (e.g., pH=3-5), the coupling is much stronger. As the potential approaches the diffusion-limited current region, a sharp drop in the local H₃O⁺ concentration (a \"jump\" in the local pH) occurs. This strong nonlinearity increases the difficulty of convergence. Under these conditions, our algorithm typically requires at least 20 iterations to achieve a self-consistent solution for the local environment and the surface kinetics across the entire polarization curve.\n2. Treatment of Hydrogen Coverage and Selection of θ\n*H\nTo clarify, the microkinetic modeling (MKM) was primarily based on energetics calculated at 0 ML hydrogen coverage. From these data, we obtained the potential- and pH-dependent H* coverage on Au(111), as shown in Figure. R3. The results reveal that the H* coverage is highly sensitive to both potential and local pH. At 0 V vs. SHE, the weak adsorption of H* on Au renders the surface essentially free of hydrogen, regardless of the pH. As the potential becomes more negative, the coverage exhibits distinct behaviors depending on the reaction condition. In acidic and neutral conditions, hydronium dissociation in the Volmer step is facile, resulting in nearly full H* coverage, although the subsequent Tafel step governs the overall kinetics. In contrast, under alkaline conditions, the alkaline Volmer step dominates with a much higher barrier, leading to negligible H* coverage. These trends are consistent with the reaction energetics and rate-determining steps discussed in the main text (Fig. 4C). To capture coverage effects more explicitly, we also calculated adsorption energetics at 7/9 ML hydrogen coverage. For intermediate coverages, we assume the adsorption energy to vary linearly with the surface H* coverage, allowing us to interpolate the kinetics parameters across the full range using the two datasets in Supplementary Tables 6 and 7.\nFigure R4. Calculated heatmap of H* coverage on Au(111) obtained from microkinetic modeling (MKM) using energetics referenced to the 0 ML state.\nFurthermore, we provide the surface Pourbaix diagram of Au(111) in Figure. R4. The diagram confirms that the equilibrium H* coverage remains very low in acidic and neutral electrolytes, while appreciable H* coverage emerges under more alkaline potentials. We note, however, that the Pourbaix diagram reflects purely thermodynamic trends and cannot capture shifts in the rate-determining step, which are responsible for the potential-dependent coverage changes seen in the MKM. Nevertheless, both approaches consistently indicate that Au(111) can host either near-zero or near-full H* coverage within the potential window studied, underscoring the importance of accounting for coverage effects in our simulations.\nFigure R5. Surface Pourbaix diagram of H* on Au(111) constructed from thermodynamic adsorption free energies.\nFollowing your valuable suggestions, we have added this crucial information to the revised SI. Please see the lines 178-182 on page S9 for “Convergence of iteration” and lines 86-88 on page S4 for “Treatment of coverage” in the revised SI and also below:\n“This MKM-transport coupling typically converges in fewer than 4 iterations in strongly acidic/alkaline conditions, but requires at least 20 iterations in dilute acid conditions (pH 3-5) due to the sharp local pH changes near the diffusion-limited region. Convergence was achieved when the relative change in local pH and H₂ concentration fell below 10-3.”\n“The microkinetic modeling utilized energetics calculated at 0 ML hydrogen coverage as reference, with linear interpolation between 0 ML and 7/9 ML datasets to account for coverage effects.”\nComment 2. The hybrid explicit/implicit model used for the GC-DFT calculation lacks supporting evidence. The authors need to discuss the origin of this model and provide a rigorous test to demonstrate its validity. This is crucial for the scientific community to trust the model as an accurate methodology for simulating solid-liquid interfaces. The current approach appears to be heavily biased by the authors' perspective.\nAuthor Reply: We appreciate the reviewer's insightful comments regarding the validation of our hybrid explicit/implicit solvation model. We agree that establishing methodological credibility is crucial for the community's trust in simulating solid-liquid interfaces.\nOur approach is grounded in established computational electrochemical practices. The implicit electrolyte description, parameterized with the bulk dielectric constant and Debye length, has been extensively benchmarked against experimental solvation free energies and validated in previous interfacial studies [J. Chem. Phys., 2019, 151]. This is self-consistently integrated with the GC-DFT formalism [J. Phys. Chem. Lett. 2016, 7, 1471-1477. J. Chem. Phys. 2017, 146, No. 114104] to explicitly control the electrode potential. This combined implicit/GC-DFT framework provides a more physicochemical and accurate description of the electrochemical interface than conventional DFT.\nTo capture specific chemical interactions, such as hydrogen bonding, proton transfer, and the solvation of intermediates, we introduced explicit water molecules and hydronium. This hybrid strategy is supported by prior rigorous work: Li et al. directly compared energy barriers between a bilayer water model and AIMD, finding differences of less than 0.1 eV across a wide potential window (JACS Au 2023, 3, 943-952), while Goddard III et. al. successfully employed 3-4 explicit waters to model proton transfer within GC-DFT (Nat. Commun. 2020, 11, 2256). These studies confirm that a minimal explicit layer can capture essential solvent effects with high fidelity at a fraction of the computational cost of fully explicit solvation.\nWe acknowledge that while ab initio molecular dynamics (AIMD) with a fully explicit electrolyte would provide a more detailed dynamic representation of the solvation environment, such an approach faces fundamental limitations for addressing our research objectives. Our continuum transport model demonstrates that transport limitations can cause variations of several orders of magnitude in the concentration of H₃O⁺ or OH⁻ ions at the electrode interface, yet standard AIMD cannot achieve precise control over such a wide range of ion concentrations, particularly the local pH that is essential for establishing quantitative free energy references. Furthermore, fully dynamic solvent configurations would substantially increase computational costs and pose convergence challenges in CI-NEB barrier calculations.\nTo address our primary research objective of understanding how mass transport effects influence reaction performance and the evolution of local environments through coupling first-principles calculations with mass transport modeling, we implemented this hybrid explicit/implicit solvation model for the electrochemical interface. This methodology provides dual benefits: the explicit solvent component captures essential hydrogen-bond interactions between intermediates and surrounding H₂O molecules, while the implicit solvation component accurately describes long-range electrostatic interactions and ensures proper decay of the electrostatic potential in the electrolyte region. This integrated strategy represents a practical and computationally efficient approach well-suited to our research goals.\nTo clarify this, we included a more detailed discussion of our model's origin, its validation against literature benchmarks, specific advantages, and limitations in the revised manuscript.\nPlease see the lines 36-41 on page S2 in the revised SI and also below:\n“The hybrid explicit/implicit solvation model builds upon established computational electrochemistry approaches [J. Chem. Phys., 2019, 151; J. Phys. Chem. Lett. 2016, 7, 1471; J. Chem. Phys. 2017, 146, No. 114104], combining implicit electrolyte screening with explicit water molecules to capture essential hydrogen bonding effects. The methodology is supported by previous benchmarks confirming its reliability for interfacial simulations [JACS Au 2023, 3, 943; Nat. Commun. 2020, 11, 2256], providing an ideal balance between accuracy and efficiency for investigating multiscale electrochemical phenomena.”\nComment 3. The statement on page 10, \" This convergence explains why polarization curves exhibit nearly identical trends at pH(local) = 5–13 when currents exceed 10−1 mA/cm2, as…\", is confusing. Based on Figure 3c, I observe a deviation from this explanation. At this current density, the local pH values for bulk pH 5 and 13 appear different (e.g., local pH 9 for bulk pH 5 and local pH 13 for bulk pH 13).\nAuthor Reply: We sincerely thank the reviewer for this insightful comment, which has provided us with a valuable opportunity to clarify the complex interplay between the local environment and reaction kinetics.\nUpon re-examination of Figure 3C, we confirm that the local pH for the bulk pH = 13 system remains stable over a wider current density range, indicating a distinct evolution path compared to the bulk pH = 5 system. We acknowledge that our original phrasing, \"convergence… with unified local pH conditions,\" did not accurately capture the details presented in Figure 3C. We apologize for this lack of precision and have revised the text accordingly.\nThis comment allows us to clarify and emphasize a more fundamental finding of our work: although the evolution paths of local pH differ under various alkaline conditions (bulk pH = 5–13), our multiscale simulations and experimental data unequivocally demonstrate that the corresponding HER polarization curves exhibit a high degree of consistency across the investigated potential window (Figure 3A). Indeed, this observation points to a mechanistic insight deeper than our initial description. For systems with bulk pH = 5-13, the local pH remains within the alkaline regime once the current density exceeds ~10⁻1 mA/cm². Across this alkaline regime, the macroscopic HER performance exhibits remarkable insensitivity to the specific value of the local pH. We infer that this is due to the shift of the rate-determining step to the reduction of water molecules under these conditions, a process whose rate is primarily governed by the applied electrode potential and shows weak dependence on the local pH. Therefore, the coincident macroscopic polarization curves result from two key factors: the convergence of the local pH to a universally alkaline regime and the intrinsic insensitivity of HER kinetics to pH under alkaline conditions.\nTo address this point directly, we have revised the relevant sentences in our revised manuscript. This revision more accurately and rigorously reflects our scientific findings. We are grateful to the reviewer for their valuable input, which has greatly enhanced the clarity and quality of our manuscript.\nPlease see the lines 224-230 on pages 10 in the revised manuscript and also below:\n“The coincident behavior of the polarization curves across a broad bulk pH range (5–13) can be attributed to the combined influence of two key conditions established in this regime: the convergence of the local pH to a universally alkaline regime once currents exceed 10⁻1 mA/cm², and the remarkable insensitivity of the HER kinetics to its precise value within that range [ChemElectroChem, 2014, 1, 1497-1507]. Consequently, the corresponding HER polarization curves exhibit a high degree of consistency across the investigated potential window.”\nComment 4. In the section titled \" Understanding the Mechanism of Polarization Curve Variation Trend under Acidic Conditions,\" the authors discuss results by setting the bulk pH to 3. It would be more impactful to show results for a strong acid, a strong base, and a neutral solution (e.g., pH 0, 14, and 7). The authors also need to rationalize why they chose this specific pH value.\nAuthor Reply: We thank the reviewer for this valuable suggestion. We agree that investigating a broader pH range will make our work more insightful and impactful. The primary goal of this section was to mechanistically deconvolute the different controlling regimes (kinetic, mixed, and mass-transport control) within a single polarization curve, which is a key advantage of our multiscale framework.\nChoosing bulk pH = 3 was deliberate and based on the following scientific rationale:\n(1) At a very strong acid (e.g., pH=0): The bulk concentration of H₃O⁺ is so high that mass transport is speedy. Consequently, within experimentally measurable current densities, the polarization curve is dominated by the kinetics of the proton-coupled electron transfer (PCET) step, and a well-defined diffusion-limited current plateau is typically not observed. We acknowledge that, in principle, a diffusion limitation might emerge at exceedingly high current densities due to the depletion of H₃O⁺ in the local environment. However, under such extreme conditions in a strong acid, the vigorous evolution of hydrogen bubbles would likely block active sites and dominate the observed current response, making the pure mass-transport limitation practically unobservable. This complex bubble-induced phenomenon falls beyond the scope of our current continuum transport model.\n(2) At a very strong base (e.g., pH=14) and at neutral conditions (e.g., pH=7): In both cases, the concentration of H₃O⁺ is exceedingly low. The proton source for the HER effectively shifts to H₂O molecules. Therefore, the reaction is inherently under PCET control (with water dissociation being a critical part of the process) across the entire potential range, and a diffusion-limited current plateau related to H₃O⁺ depletion is absent. As shown in Fig. 3B, the experimental results indeed demonstrate that the polarization curves at pH=7 and pH=14 exhibit very similar trends, dominated by kinetic limitations [Nat. Chem., 2013, 5, 300-306].\n(3) At pH = 3 (a dilute acid): This is an ideal situation where all three controlling regimes—kinetic, mixed, and distinct mass-transport control (evidenced by a clear current plateau due to H₃O⁺ depletion)—are prominently visible within a standard experimental window. This makes it an ideal and non-arbitrary case study for demonstrating how our framework can dissect the origins of each segment of the polarization curve. The existence of distinct controlling regimes (PCET → mixed → mass-transport→ mixed→ PCET) in dilute acids is well-documented experimentally. Our choice is strongly supported by the work of [J. Phys. Chem. C, 2020, 124, 3988−4000], which systematically analyzed HER polarization curves on Au and Pt across a dilute acid range (pH 2.5–3.9). That study explicitly identified, through experimental data and analytical modeling, the characteristic \"two-step\" behavior of the polarization curve where the controlling factor shifts with potential. Our multiscale modeling builds upon this foundational experimental observation by providing an atomistic-scale resolution and a theoretical explanation for the underlying mechanisms governing these transitions. We will cite this work appropriately in the revised manuscript to ground our rationale firmly in established experimental observations.\nPlease see the lines 321-334 on page 14 in the revised manuscript and also below:\n“Notably, the choice of pH = 3 for our detailed analysis of polarization behavior under acidic conditions was guided by fundamental electrochemical considerations across the pH range. In strongly acidic media, the high bulk concentration of H O⁺ leads to rapid mass transport, resulting in polarization curves dominated by proton-coupled electron transfer kinetics without a well-defined diffusion-limited current plateau. In alkaline and neutral environments, the limited availability of H O⁺ shifts the proton source to water molecules, maintaining the reaction under kinetic control throughout the potential range. In contrast, dilute acid conditions at pH = 3 provide a unique window where different controlling regimes are prominently exhibited within experimentally accessible current densities. This behavior is well-documented in systematic studies of HER polarization curves [J. Phys. Chem. C, 2020, 124, 3988-4000]. Our multiscale framework builds upon these experimental observations by providing atomistic-scale insights into the mechanisms governing these transitions, making pH = 3 an ideal case for demonstrating how our model deconvolutes the origins of different polarization curve segments.”\nComment 5. Although the authors discuss their methodology, they do not provide any input/output files needed to reproduce their results. To ensure reproducibility and transparency, they must share the raw input and results data for the GC-DFT, MKM, and transport modeling on a public platform like a Git Hub page. This is essential and should be shared during the review process.\nAuthor Reply: We thank the reviewer for this essential suggestion regarding research reproducibility. In response, we have created a comprehensive GitHub repository containing all input and output files necessary to reproduce our results. The repository includes:\n• Coordinates and calculated results for all adsorbates and transition states of GC-DFT\n• Microkinetic model source code and parameters\n• Continuum transport model implementation\n• Detailed documentation and setup instructions\nThe repository is now publicly available at: https://github.com/WestlakeUniversity-WangLab/Multiscale_HER.\nWe have added the repository link to the Data Availability section in the revised manuscript and will maintain this resource for community access.\nReviewer #3 (Remarks to the Author):\nThe paper “Toward Rational Understanding of the Hydrogen Evolution Polarization Curves Through Integrating Ab Initio Simulations and Mass Transport Effects” presents a valuable modeling study aimed at understanding polarization curves for the hydrogen evolution reaction (HER) on Au(111) and Pt(111) surfaces. The authors employ multiscale simulations that integrate grand-canonical density functional theory (GC-DFT) calculations, a microkinetic model, and a continuum mass transport model.\nThe central contribution of the paper is the demonstration that mass transport phenomena significantly influence polarization behavior and catalyst performance, and that their explicit inclusion in simulations is essential to accurately reproduce experimental polarization curves.\nOverall, the work addresses a timely and important topic in the theoretical description of electrochemical processes. The manuscript contains sufficient novelty to merit consideration for publication in Nature Communications. However, several issues should be addressed prior to acceptance:\nAuthor Reply: We thank the reviewer for the professional vision and constructive suggestions, which greatly enhance the quality of our work. Following your valuable suggestions, we revised our manuscript, which we sincerely hope will address your concerns. Thank you for your professional input, which has helped us shape the manuscript into a more comprehensive and impactful study.\nComment 1. In the Abstract and Introduction, the authors state that results are presented for both Au(111) and Pt(111). However, the Results section primarily discusses Au(111). I recommend either including Pt(111) data in the main text or removing its mention from the Abstract.\nAuthor Reply: We thank the reviewer for this pertinent observation. In our initial submission, we referenced Pt(111) primarily to demonstrate the robustness and reproducibility of our multiscale framework using previously published GC-DFT data [J Am Chem Soc, 2020, 142, 4985-4989]. However, as the reviewer correctly notes, the detailed results and discussion in the manuscript focus predominantly on Au(111), for which we conducted original GC-DFT calculations, enabling deeper mechanistic analysis.\nTo maintain consistency between the stated scope and the actual content, and to avoid potential confusion, we have removed the mention of Pt(111) from the Abstract and Introduction in the revised manuscript. This allows us to concentrate fully on presenting a coherent and comprehensive analysis of the Au(111) system, for which we can provide complete methodological transparency and in-depth discussion of the underlying reaction energetics and kinetics.\nComment 2. Many sentences in the Abstract and Introduction are repetitive. Please revise to improve conciseness and avoid redundancy.\nAuthor Reply: We thank the reviewer for identifying repetitiveness in our Abstract and Introduction. We have comprehensively revised both sections to enhance conciseness and eliminate redundancy.\nThe Abstract has been streamlined by removing repetitive methodological descriptions and consolidating overlapping statements about mass transport effects. The Introduction has been restructured to eliminate redundant discussions of catalyst design challenges. The literature review now presents a concise analysis of methodological evolution rather than listing previous works. Our framework description and achievements have been integrated into a cohesive narrative that emphasizes scientific discovery over repetitive claims. These revisions have significantly improved the clarity and flow of both sections while maintaining all essential scientific content.\nComment 3. The reference energy of the products used in the free energy diagrams (Fig. 2a and Fig. 2b) is not clearly defined. The method of determining these reference energies should be explicitly stated.\nAuthor Reply: We thank the reviewer for this insightful comment and apologize for the lack of clarity regarding the reference energies in our free energy diagrams (Figs. 2a and 2b).\nThe energy profiles were constructed along the reaction coordinate, with the energy zero (0 eV) defined as the initial state of the first elementary step (the Volmer step) for each pathway. Specifically:\n• For the acidic pathway, the reference is: * + H₃O⁺ + e⁻\n• For the alkaline pathway, the reference is: * + H₂O + e⁻\nAll subsequent intermediates and transition states are calculated relative to this common initial state. For example, in the Volmer-Heyrovsky pathway, the energy of the H* intermediate generated in the Volmer step serves as the reference for calculating the energy change of the subsequent Heyrovsky step. This approach maintains thermodynamic consistency across the entire reaction landscape.\nTo systematically evaluate the effect of hydrogen coverage on the reaction energetics, we constructed separate energy profiles for two distinct initial surface environments:\n• θ *H = 0 ML: clean surface without pre-adsorbed H\n• θ *H = 7/9 ML: surface with seven pre-adsorbed H atoms in a unit cell\nThe corresponding atomic structures are provided in Supplementary Figures 5 and 6.\nWe will revise the manuscript to explicitly state this reference energy definition and the stepwise construction methodology for the free energy diagrams.\nPlease see the lines 69-78 on pages S3-S4 in the revised SI and also below:\n“Free Energy Diagrams\nThe free energy diagrams (Figs. 2A, B) use the initial Volmer step state as energy zero (0 eV):\n• For the acidic pathway, the reference is: * + H₃O⁺ + e⁻\n• For the alkaline pathway, the reference is: * + H₂O + e⁻\nAll subsequent states are calculated relative to this reference. To evaluate coverage effects, we compared two surface environments:\n• θ *H = 0 ML: clean surface without pre-adsorbed H\n• θ *H = 7/9 ML surface with seven pre-adsorbed H atoms in a unit cell\nCorresponding structures are shown in Supplementary Figures 5-6”\nComment 4. The definition and calculation procedure of the Volmer step ratio (%) reported in Fig. 4b are unclear and require clarification.\nAuthor Reply: We thank the reviewer for this insightful comment, which allows us to clarify the definition and significance of the Volmer step ratio. We apologize for the unclear description and have clarified this in our revised manuscript.\nMore specifically, the Volmer step ratio is defined as the percentage contribution of the acidic Volmer step to the total production rate of the H* intermediate. It specifically distinguishes the proton source for the H* formation:\n• Acidic Volmer: H₃O⁺ + * + e⁻ → H* + H₂O\n• Alkaline Volmer: H₂O + * + e⁻ → H* + OH⁻\nThe ratio is calculated from the steady-state solution of our microkinetic model (MKM), which is fed the local environment from the continuum transport model. At a given potential, the MKM solves a system of ODEs for surface coverages. The reaction rates of the acidic (r ) and alkaline (r ) Volmer steps are direct outputs of this steady-state solution. The ratio is then computed as: r / (r + r )] × 100%.\nThis metric quantitatively reveals the dominant proton source (H₃O⁺ vs. H₂O) for the HER under different potentials. As shown in Fig. 4b, the decrease of this ratio with increasing negative potential signifies the depletion of H₃O⁺ near the electrode due to mass transport limitations. When the ratio approaches zero, it indicates that H* is exclusively generated from water dissociation, and the reaction re-enters a regime dominated purely by the kinetics of the Proton-Coupled Electron Transfer (PCET) step.\nWe have included a clear definition of the Volmer step ratio and a concise description of its calculation in our revised manuscript.\nPlease see the lines 118-127 on page S6-7 in the revised SI and also below:\n“Definition and Calculation of Volmer Step Ratio\nThe Volmer step ratio quantifies the relative contribution of the acidic versus alkaline proton transfer pathways to the formation of the H* intermediate. It is defined as the percentage of the total H* production rate that proceeds through the acidic Volmer mechanism:\nVolmer step ratio = r / (r + r) ×100%\nThe individual reaction rates are obtained from the solution of the microkinetic model, which is self-consistently coupled with the continuum transport calculations. The potential dependence of this ratio reflects the shifting availability of H O⁺ versus H O as proton donors under different interfacial environments.”\nComment 5. The manuscript introduces the concept of “locality” of the pH. How is this defined? If it corresponds to the concentration of H₃O⁺ within the boundary layer, please provide the boundary layer thickness, at least for the case of an RDE rotation speed of 1600 rpm.\nAuthor Reply: We thank the reviewer for this insightful suggestion, which allows us to clarify the connection between our model and the experimental conditions.\nThe reviewer's understanding is correct. In our study, the local pH is defined as the pH value at the cathode surface. It is directly calculated from the local concentration of hydronium ions (c ) at that plane using the standard formula:\npH = −log (c )\n(local) 10 H3O+,local\nwhere c has units of mol/L. As schematically illustrated in Supplementary Figure 1, c denotes the concentration of H O⁺ ions at the cathode surface. This c is a key output of our coupled multiscale model. It results from the balance between the consumption of H O⁺ by the surface reaction and its supply by diffusion, migration, and convection from the bulk solution. The characteristic length scale for this transport process is indeed the diffusion boundary layer thickness (δ).\nRegarding the specific value requested by the reviewer, the boundary layer thickness used in our simulations is δ = 4 × 10⁻4 m (400 μm). This parameter is stated in the Methods section of the Supplementary Information. In summary, the local pH is a direct consequence of the reaction-transport balance across this defined boundary layer of thickness δ.\nWe have included a specific definition of local pH for clarity in our revised manuscript and the thickness of boundary layer (δ) was listed in Supplementary Table 2.\nPlease see the lines 83-86 on page S4 in the revised SI and also below:\n“The local pH is defined as the local pH at the cathode surface, which is directly calculated from the local concentration of hydronium ions (c ) at that plane, as schematically illustrated in Supplementary Figure 1.”\nComment 6. The flow of argumentation in the Results section is at times difficult to follow. Improving the clarity and organization of this section would strengthen the paper.\nAuthor Reply: We sincerely thank the reviewer for this crucial comment regarding the organization of our Results section. We have completely restructured this section to enhance logical flow and clarity. The revised Results section now features a clearer narrative progression through the following reorganized subsections:\nWe clarify the Results section into 6 parts in a logical way:\n1. Multiscale Simulation Framework: We begin by clearly presenting our computational methodology, establishing the foundation for all subsequent analysis.\n2. Reaction Mechanisms and Free Energy Diagram of HER on the Au(111) Surface: With the framework established, we elucidates the atomic-scale fundamentals, presenting the reaction pathways and free energy diagrams.\n3. Multiscale Simulations of Polarization Curves of HER on the Au(111) Surface: This section validates our multiscale framework by showing its predictive power by directly comparing simulated and experimental polarization curves across the entire pH range, bridging from microscopic GC-DFT calculations to macroscopic electrochemical performance.\n4. Understanding the Effects of pH on Polarization Curves: Building on the validated model, this section establishes the significant deviation of the local environment from bulk conditions. We identify distinct potential-dependent regimes of local pH and elucidate their crucial influence on the macroscopic HER polarization behavior.\n5. Understanding Mechanisms behind the Variation Trend of the Polarization Curves under Acidic Conditions: By providing deep mechanistic insights into the polarization behavior under acidic conditions, theoretically classifying it into five distinct regimes governed by different factors.\n6. Fine-tuning and Optimizations of Mass Transport: This section validates the five-regime classification by systematically varying catalyst types and rotation rates to represent modifications in PCET and mass transport effects. We further elucidate how rotation rate adjustments influence the local environment and define the operational limits of such mass transport control.\nBeyond the structural reorganization, we have refined the textual content of each parts to enhance readability and logical flow and the key results discussions have been refined to strengthen their connections with corresponding data analyses. This redesigned structure ensures a coherent progression from methodology establishment to model validation, mechanistic exploration, and practical application.\nWe believe these revisions have substantially improved the manuscript's clarity and organization, and we sincerely thank the reviewer for this invaluable guidance in facilitating these enhancements.\nComment 7. The expression for the time derivative of H* coverage (Eq. 15) is formally incorrect. The two terms associated with the Tafel step should be multiplied by a factor of 2 (see Eq. 115 in [Int. J. Hydrogen Energy, 2019, DOI: 10.1016/j.ijhydene.2019.05.183]). Please verify whether this is a typographical error or requires correction.\nAuthor Reply: We sincerely thank the reviewer for your meticulous reading and for pointing out this error in Equation 15. The reviewer is correct. The two terms associated with the Tafel step should indeed be multiplied by a factor of 2, as this step (H* + H* → H₂) consumes two H* species, a standard requirement correctly referenced from the literature [Int. J. Hydrogen Energy, 2019, DOI: 10.1016/j.ijhydene.2019.05.183]. This was a typographical error in the manuscript text.\nThe corrected Equation (15) is as follows:\nWhere the two terms associated with the Tafel step now correctly account for the consumption of two H* species in the Tafel step.\nWe have corrected Equation 15 in the revised manuscript to accurately represent the rate of change of H* coverage. Furthermore, to ensure complete transparency and reproducibility, we double-checked all equations in the Supplementary Information and updated them as needed to guarantee consistency. We confirm that this was a typographical error, and the correct formulation (accounting for the stoichiometric factor of 2) was implemented in all our actual microkinetic model calculations. Therefore, this typo does not affect any of the simulated results or conclusions presented in the study.\nComment 8. The manuscript reports Eq. 21 as j(H₂) = 2·F·r(H₂), with r(H₂) = r₂ + r₄ + r₅. However, r(H₂) is originally defined as r(H₂) = r₁ + r₂ + r₃ + r₄, which can be reformulated as the authors present only under the steady-state assumption (dθ(H*)/dt = r₁ – r₂ + r₃ – r₄ – 2·r₅ = 0). The authors should clarify whether the steady-state approximation has been explicitly applied in the microkinetic model.\nAuthor Reply: We sincerely thank the reviewer for this professional and insightful comment, which allows us to clarify a key point in our derivation. The reviewer is correct that the equivalence between the expressions relies on the steady-state approximation (SSA), which is indeed a fundamental component of our model for simulating steady-state polarization curves.\nWe apologize for the lack of clarity in the manuscript. The most rigorous way to derive the current density is from the perspective of electron conservation.\nElectron Conservation and Current Definition.\nThe rigorous starting point is electron conservation: the total Faradaic current equals the rate of electron consumption by all proton-coupled electron transfer (PCET) steps. In our network, these are the Volmer (r₁, r₂) and Heyrovsky (r₃, r₄), each involving one electron. Thus, the total current density can be written as,\nj(H₂) = F(r₁ + r₂ + r₃ + r₄)\nFor the H* intermediate, the steady-state approximation condition is\n(dθ(H*)/dt = r₁ – r₂ + r₃ – r₄ – 2r₅ = 0\nwhich gives\nr₁ + r₃ = r₂ + r₄ + 2r₅\nSubstituting this relation into the current expression yields\nr₁ + r₂ + r₃ + r₄ = 2(r₂ + r₄ + r₅)\nand therefore,\nj(H₂) = F(r₁ + r₂ + r₃ + r₄) = 2F(r₂ + r₄ + r₅) = 2Fr(H₂)\nDefinition of r(H₂):\nPhysically, only the Heyrovsky (r2, r4) and Tafel (r5) steps directly generate molecular hydrogen. Therefore, the expression of r(H₂) = r₂ + r₄ + r₅, which was used in the manuscript, is not the fundamental definition of the H₂ production rate, but an equivalent rate expression that is rigorously valid under the SSA. It is defined such that when multiplied by 2F, it yields the correct total current derived from electron conservation.\nWe have revised the manuscript to explicitly state the use of the steady-state approximation in the microkinetic model and include this clear derivation to explain the origin of Eq. 21 and the definition of r(H₂).\nPlease see the lines 105-115 on page S5-S6 in the revised SI and also below:\n“the formation rate of H was calculated as,\nr = r + r + r\nwhere this expression represents the effective production rate under the steady-state approximation, as rigorously derived from electron conservation and H* coverage stationarity. The corresponding consumption/production rates of other species are given by:\nr = −r −r\nr OH* = r ? + r @\nr = r + r −r −r\nr K+ = r ClO 4 * = 0\nand the partial current density of H was given by,\nj = 2Fr\nwhere the factor of 2 accounts for the two electrons required to produce each H₂ molecule, consistent with the electron balance under steady-state conditions. The flux continuity at the electrode surface will be described by,\nsj = −N / nF”"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have reasonably addressed all my previous comments and concerns. Therefore, I can recommend this work for publication. Meanwhile, I suggest that the authors enrich the descriptions of the advantages and challenges of the reported framework in tackling more complex electrochemical reactions in the discussion part. Further review is not needed.", "reviewer_2": "The authors well addressed my comments. I'm happy to see this paper in Nature Communcations online soon.", "reviewer_3": "The authors have extensively and appropriately addressed the points raised in the original reports. The paper can now be accepted for publication."}, "author_response": ""}]

Decision letter