Modeling Study of pH Distribution and Non-Equilibrium State of Water in Hydrogen Evolution Reaction

pH distribution and non-equilibrium state of water in hydrogen evolution reaction are studied using continuum models. In the ﬁrst model, we analyze the pH distribution in a rotating ring disk electrode system, where the hydrogen evolution occurs on the disk electrode. The model predicts a pH distribution comparable to the experimental data and the nonequilibrium state of water (cH ∗ cOH > 1.0 × 10 − 14 ) in a small portion of the diffusion layer (ca. 5 μ m in thickness) adjacent to the bulk electrolyte under forced convection. The second model explores the pH distribution on an electrode with nanovoids in hydrogen evolution reaction in an acidic media. The value of cH ∗ cOH shifts signiﬁcantly when close to the electrode surface, e.g., ≤ 2 μ m, indicating pH is not viable to assess its impact on an electrochemical reaction involving hydroxide ions. Modeling results also prove that, for an electrode with nanovoids, the concentration gradient of hydroxide between the plain ﬁeld and the bottom of the nanovoid is minimal. Therefore it should not be the root cause for the differential kinetics of metal electrodeposition inside/outside the nanovoids. © The This an open access the terms of the Creative Commons unrestricted [DOI: 10.1149/2.0312001JES] JES Issue on Mathematical Modeling of Electrochemical Systems at Multiple Scales in Honor of Richard Alkire.

The pH value plays a key role in many aqueous chemical/electrochemical systems, including biological tissues, chemical synthesis, electrocatalysis, and electrochemical deposition. [1][2][3][4] Reactions involving protons or hydroxides are affected by pH in both kinetics and thermodynamics, e.g., redox potential predicted by Nernst equation in electrochemical reaction. 5 A series of key electrochemical reactions for energy storage and biological systems, including hydrogen evolution reaction (HER), Oxygen Evolution Reaction (OER), Oxygen Reduction Reaction (ORR), and Hydrogen Oxidation Reaction (HOR), have been studied intensely for last decades. [6][7][8][9][10] Production or consumption of proton/hydroxide in these reactions can unambiguously alternate both the interfacial and the surrounding pH values, and ultimately shift thermodynamics and kinetics thereof. 11 Unlike the bulk aqueous electrolyte, where the pH could be readily determined by potentiometric glass electrodes etc., interfacial pH measurements typically require in situ techniques, e.g., Spectroscopy, Rotating Ring Disk Electrode (RRDE), and Scanning Electrochemical Microscopy (SECM). [12][13][14] Recently, Zimer et al. reported a versatile RRDE method to measure interfacial pH, in which the ring electrode was modified with pH-sensitive IrO x to detect the pH variation induced by the HER on the disk electrode. 15 In the RRDE system, the product of electrochemical reactions from disk can instantly transport to the ring in a fraction of one second to tens of milliseconds, depending on the rotating rate. In the absence of any homogeneous chemical reactions in the electrolyte, the flux of product to the ring with respect to that toward the disk, i.e., collection efficiency, is determined by the geometry of the RRDE assembly. 5 However, if homogeneous chemical reactions are involved, the product diffusing away from the disk electrode surface reacts fast and easily, which can distort the quantitative measurements on the ring electrode. 7,8 On the other hand, pH distribution inside nanovoids is of interest since chemical/electrochemical reactions inside these semi-closed space are typically different to those at the open surface exposed to the solution. 16,17 One example in electrochemical systems is the patterned electrode, commonly used in medical devices, sensors, and integrated chip fabrication. [18][19][20][21][22] Recent work showed the possible impact of pH on the damascene cobalt electrochemical plating in nano-trenches. 23 The pH, or the concentration of hydroxide, would affect the cobalt plating both thermodynamically and kinetically. 24 However, the actual mechanism of cobalt bottom-up filling in nanovoids is still under debate, 25 and, whether the concentration gradient of the hydroxide between the plain electrode field and the nanovoids is accountable for the differential kinetics should be examined.
Herein, we developed two continuum models to study the pH and non-equilibrium state of water in electrochemical systems: RRDE and electrode with nanovoids. The first model, diffusion-convectionreaction model, examines the RRDE method for in situ pH measurements, and reports the significant pH variation between the ring electrode and disk electrode, and also the break of water equilibrium in a small portion of diffusion layer adjacent to the bulk solution. Following that, our second model reports significant variation from water equilibrium when close to the electrode surface, e.g. ≤ 2 μm, in the case of hydrogen evolution in acidic media. Thus, instead of pH, the concentration of hydroxide can be utilized to assess its impact on an electrochemical reaction involving hydroxide ions. Further evidence shows the concentration of hydroxide is not responsible for the differential kinetics between plain electrode field and bottom of the nanovoids in Co bottom-up filling.

Model Development
Model A: rotating ring disk electrode A two-dimensional axis symmetrical diffusion-convectionreaction model (Scheme 1) presented in this work was performed Scheme 1. Schematic diagram of the two-dimensional axis-symmetrical model used in the simulation: RRDE.

Reaction
Reaction formula Description within the COMSOL Multiphysics (version 5.3) using the same conditions as the electrochemical experiments described by Zimer et al. 15 Briefly, the length and height of the domain in the model are set at two times of the outer radius of the ring electrode, R 3 = 0.362 cm, and three times of the diffusion boundary layer thickness, δ. δ is a function of the rotation rate and can be approximated by 5 where ω is the rotation rate (rad s −1 ), equal to 2π f 60 ( f , the rotation rate in the unit of rpm), D is the diffusion coefficient and ν is kinematic viscosity of the electrolyte.
The governing equation for the diffusion-convection-reaction model is 26 Where D is diffusion coefficient; c is concentration (cH and cOH represent the concentration of proton and hydroxide, respectively); z is vertical axis; r is radial axis; u r u z are fluid velocity close to the RRDE according to Levich formulism: 5 u r = 0.51023 · (ω 3 /ν) 1/2 rz [3] u z = −0.51023 · (ω 3 /ν) 1/2 z 2 [4] R H 2 O is the reaction rate of homogeneous water dissociation/association in the bulk electrolyte, as specified in R7, Table I, and is a function of the actual product of cH and cOH: 27 where k f is the rate of the forward reaction of R7. Boundary 1 is the disk electrode, where the hydrogen evolution reaction occurs through (R2) and produces hydroxide. Note that (R1) is not considered in the model here for simplification, since the reactant in hydrogen evolution in neutral pH is water (R2) as reported by Shinagawa et al. 28 The boundary condition is where the flux of hydroxide is assumed uniformly distributed on the disk electrode, considering as a kinetically controlled electrochemical reaction. As shown in Scheme 1, Boundary 2 and 4 are insulators, and no flux occurs on these two boundaries; Boundary 3 is the ring electrode and designed to measure the open circuit potential; therefore no flux occurs; Boundary 5 is the outflow, i.e. n · D∇c = 0; Boundary 6 is the bulk electrolyte, where cH = 1 × 10 −6 M; cOH = 1 × 10 −8 M; Boundary 7 is the symmetrical axis. Initial conditions in the domain are set at cH = 1 × 10 −6 M; cOH = 1 × 10 −8 M Model B: electrode with nanovoids A two-dimensional axis symmetrical diffusion-reaction model was built, as shown in Panel A Scheme 2. The governing equation for the diffusion model is 26 All the parameters in Equation 7 have the same meaning as those in Equation 2. The dissociation/association reaction of water is the same as (R7) in the first model; thus, the same kinetics, Equation 5, applies to R H 2 O in the Equation 7.
Initial conditions in the domain are set at cH = 1×10 −3 M; cOH = 1 × 10 −11 M as those specified in the experimental work. 23 Boundary 1 is the working electrode (R 1 = 1 μm); No flux occurs on Boundary 2; Boundary 3 was the edge of the diffusion layer, assuming the thickness of the diffusion layer, δ, is 150 μm; Boundary 4 is the symmetrical axis. For the electrode with nanovoids (Panel B, Scheme 2), all boundaries were the same as those for the plain electrode, except a nanovoid is present (d = 100 nm, R 2 = 10 nm). The bottom of the nanovoid is labeled as "bottom" and the plain part of the electrode is labeled as "field". Line A, perpendicular to the bottom of the void, and Line B, perpendicular to the field of the plain electrode are highlighted in Panel B, Scheme 2, to study the concentration distribution inside and outside the nanovoid. In both cases, the hydrogen evolution reaction is assumed to occur through proton reduction (R1) and the concentration of protons on the electrode is assumed to be 0. 29 All parameters described in both model A and B above are listed in Table II.

Results and Discussion
Rotating ring disk electrode.-Multiple current densities digitized from Chronoamperometry with a rotation rate at 1800 rpm, 15   i. The pH drops quickly along the RRDE surface and is much lower on the ring electrode (r = 3.61mm) compared to that on the disk electrode. It clearly shows the inaccuracy to assume the interfacial pH on the disk electrode from the measurements on the ring due to the fast water association reaction (R7) in the diffusion layer which is as thin as tens of micrometers. ii. While the pH values calculated with 1 * D OH are much more alkaline than experimental values, the ones calculated with 5 * D OH compare well with the experimental data (Green Squares) in both panels. It is in good agreement with the recent report by Grozovski et al., 11 in which hydroxide is found to diffuse much faster in the highly non-equilibrium state through a "directed Grotthuss mechanism", i.e., a shuttling of O-H bond responsible for the large value of D OH . 16,30,31 iii. A significant discrepancy between experimental data and modeling results appears in Panel C, Figure 1, where the higher current is flowing on the disk electrode. This phenomenon could be explained as following: Open circuit potential, OCP, of the IrO x ring electrode, converted into pH through a standard curve, is determined by both the electrostatic potential and the pH of the electrolyte in the actual experiment. Large magnitude cathodic current flowing in the electrolyte, as i c in Panel C (Figure 1), may yield a significant negative electrostatic potential (smaller OCP) monitored by the ring electrode, and therefore, reporting a higher pH value. [32][33][34] Also, the IrO x ring electrode shows memory effect in the actual experiment, e.g., OCP of IrO x shifted positively after experiencing much more alkaline media as described in the experimental work. The phenomena are especially pronounced with the Chronoamperometry experiment, where the initial current producing hydroxide is as high as −5 mA cm −2 before it achieves the steady-state, i c = −1.35 mA cm −2 . As shown in the experimental work, 15 the offset of OCP is c.a. 9.4 mV before and after the single measurement, corresponding to a shift of 0.4 in pH, which is not captured by the model.
For electrochemical redox species not involving bulk reactions, the thickness of the diffusion layer on the rotating disk electrode could be calculated by Equation 1. δ OH is calculated to be 35 μm for 5 * D OH . However, the model shows that the diffusion layer thickness increases from 25 to 40 μm as more hydroxide is produced (see Figure 2, Panel A). Unlike the linear concentration profile of electrochemical reactant/product in the absence of bulk chemical reactions in the diffusion layer, 5 the concentration gradient of both the proton and hydroxide is steeper in the part of the diffusion layer adjacent to the bulk electrolyte (Panel A and B, Figure 2). A similar concentration gradient profile was reported by Grozoyski et al. 11 in which a constraint of cH * cOH = 1 × 10 −14 M 2 was employed to the whole domain.
Because of the fast kinetics of acid-base reactions, water is typically assumed to be in equilibrium, e.g., cH * cOH = 1 × 10 −14 . However, it is found an unordinary region where cH * cOH> 1 × 10 −14 could be obtained in the part of diffusion layer adjacent to the bulk electrolyte as shown in Panel C and Panel D,   hydroxide decreases cH * cOH. At 5 * D OH and cH * cOH = 1 × 10 −14 is well obeyed with smaller current density, I a , in the whole domain.
The unique phenomenon could be attributed to the enhanced mass transport by the steep concentration gradient of both ionic species, which only achieves steady-state with faster reaction rate of acid-base bulk reactions. As a matter of fact, the term determining the concentration profile of ionic species in mass balance equation (Equation 2) is mass transport and the overall reaction rate (Equation 5). A careful inspection of the Equation 5 will show that the term in parenthesis, Kw-cH * cOH, will generate a small value on the order of 10 −14 M 2 . Thus the mass transport of protons and hydroxide becomes comparable under the steep concentration gradient. The interpretation above is proven by the data in Figure 2, i.e. for same current density, i a , a larger diffusion coefficient (5 * D OH vs 1 * D OH ), produces a smaller concentration gradient (black curves in Panel A and Panel B), thus cH * cOH is maintained at 1 × 10 −14 at in the whole domain (black curve in Panel C) while cH * cOH increases above 1 × 10 −14 at 1 * D OH (black curve in Panel D). To further confirm our explanation for the phenomenon, the rate constant of water association reactions, k f , is tuned up to 3 orders of magnitude higher, and a constant cH * cOH, i.e. 1 × 10 −14 M 2 is well achieved in all scenarios (Data not shown). This observation affords strong evidence that the non-equilibrium state of water in the diffusion layer adjacent to the bulk electrolyte originates from the competition between mass transport and homogeneous reaction rate of both protons and hydroxides. The pH profiles along the lines normal to the center of the disk electrode (solid lines) and the inner edge of the ring electrode (dashed lines) are depicted in Figure 3. A higher flux of hydroxide yields more alkaline environment, e.g., pH = 9.4 for I a , pH = 9.8 for I b , and pH = 10.3 for I c on the disk electrode surface, corresponding to pH = 7.2, 8.5, and 9.4 on the ring electrode surface. The difference in pH between that along the center line and that along inner ring edge line is most significant at the smallest cathodic current density, I a (see black lines), while they are closer at the largest current density, I c (blue lines). This phenomenon shows the limited buffering capacity of water without any additional weak acid/conjugated base.
Assuming same flux of hydroxide to the working electrode, I a , the pH on the center of the disk electrode, point A on the boundary 2 (0.25 cm from the center of the disk), and the ring electrode is examined at different rotation rates as shown in Figure 4. The pH on the center of the disk electrode (black connected triangle) and point A (black connected sphere) only decrease slightly by increasing rotation rate from 300 rpm to 1800 rpm, while the pH on the ring decreases sharply from 8.5 at 300 rpm to 7.5 at 900 rpm and gradually to 7.1 at 1800 rpm. Even though the mass transport of the produced hydroxide from the disk electrode to the ring is enhanced with a higher rotation rate, the electrolyte becomes less alkaline in this case primarily due to a larger flux of hydroxide/proton in the thinner diffusion layer at higher rotation rates, as clearly shown in the hydroxide concentration profile (Panel B, Figure 4). pH evolution in nano-features.-pH evolution in nano-features has been of interest in both research and applications due to the special chemical environment originated from the limited mass transport to the confined spaces. In the case here, we modeled pH distribution inside nanovoids of metal electrodes by assuming the limiting current of HER controlled by diffusion, i.e., the concentration of protons on the electrode surface is zero. As indicated by Grozovski  occurs primarily through proton reduction (R1) in a mildly acidic electrolyte. 11 The model starts with the plain electrode, as shown in panel A, Scheme 2 (diagram is not proportional to the actual dimension for clarity). Shown in panel A, Figure 5, is the water equilibrium along the symmetrical axis calculated with 1 * D OH (black line) and 5 * D OH (red line). Significant deviation of cH * cOH from 1 × 10 −14 is found in the region below 2 μm for both curves, which is close to 0 at the electrode surface. The non-equilibrium state of water in model B could be explained by the fast drain of protons from the diffusion layer adjacent to the electrode. As a matter of fact, a calculation shows the flux of proton to the plain electrode is 2.0 × 10 −16 mol s −1 , while the generation of protons and hydroxide, assuming only dissociation reaction occurs (R H 2 O =k f * Kw), is only 8.8 × 10 −18 mol s −1 in the domain of 2 μm height (volume = 2μm × π(R 1 ) 2 ).
Concentrations of hydroxide calculated with 1 * D OH (black line) and 5 * D OH (red line) are elucidated in panel B, Figure 5. In the region close to the electrode surface, concentration of hydroxide calculated with larger diffusion coefficient, 5 * D OH , is more diluted compared these calculated with 1 * D OH ; the hydroxide concentration profiles obtained after ca. 1 μm are essentially the same for these two parameters. pH, defined as − log(cH ), is plotted in Panel C, Figure 5. The pH profile calculated with 1 * D OH (black dashed line) and 5 * D OH (red dotted line) are virtually the same, implying one needs to consider the concentration of hydroxide, not pH, in the case of chemical/electrochemical reactions involving hydroxide, as illustrated by Co deposition reported elsewhere. 24 The electrode with nanovoid, in the dimension of depth, d = 100nm, and radius, R 2 = 10nm, is depicted in Panel B, Scheme 2. Figure 6 shows the cH * cOH, cOH, and pH along the Line A and B calculated with 1 * D OH (left column) and 5 * D OH (right column). Similar to the plane electrode, cH * cOH, is significantly away from the equilibrium value, 1 × 10 −14 . Specifically, the value achieves close to 0 in the domain of the nanovoid (blue circles in Panel A, Figure 6). The hydroxide concentration profiles are the same for Line A and B beyond the field of the electrode. The concentration of hydroxide is 1.  Figure 6). The pH inside the nanovoid goes up to 12 at the bottom of the nanovoid compared to 8 on the field (Panel C, Figure 6). A plateau of pH inside the nanovoid, ca. 11.5, is observed as highlighted in the expanded region from the blue circle. A similar phenomenon is found with 5 * D OH calculation. The major difference is that the concentration of hydroxide is lower when calculated with 5 * D OH (Panel E, Figure 6), compared to that calculated with 1 * D OH as shown in Panel B, Figure 6. The concentration of hydroxide at the bottom of the nanovoid and the field are even closer in the case of 5 * D OH , 6.6 × 10 −9 vs. 6.3 × 10 −9 M.
It should be stressed that these findings contradict the proposed mechanism for damascene cobalt electrochemical deposition, where pH is claimed to be the major driving force for the pronounced kinetics at the bottom of the nanovoid. 23 Firstly, it is inappropriate to utilize pH in the highly non-equilibrium electrolyte to demonstrate an electrochemical reaction involving hydroxide where cH * cOH is a spatially resolved factor. Electrochemical plating of cobalt is one of the examples where the formation of Co(OH) is the rate-determining step as reported elsewhere. 24,35,36 Secondly, the difference in hydroxide concentration between the field and the bottom of the nanovoid is negligible to account for the high Co growth kinetics in the void. Meanwhile, another work by Wu et al. reported a different mechanism in virtually the same electrochemical system: 25 the electrochemical reduction of the suppressor by hydrogen radical in the nanovoids is much faster than that on the field during the Co electrochemical plating process, which ultimately leads to lower surface concentration of suppressor in the nanovoid and thus faster Co reduction.

Conclusions
The major findings emerging from the simulations presented in this work relate to the interfacial pH distribution and non-equilibrium state of water in electrocatalysis involving RRDE and electrode with nanovoids via continuum models.
For RRDE, the pH distribution in the water-splitting reaction is captured by the electrochemical model and the pH value on the ring   compares well to the experimental value. The electrolyte adjacent to the ring electrode is much less alkaline from that on the disk electrode in the presence of ultrafast water association reaction. Non-Equilibrium State of water, cH * cOH> 1 × 10 −14 , is found in the sub-region of the diffusion layer which is adjacent to the bulk electrolyte, and also, the extent of the deviation from the regular water equilibrium, 1 × 10 −14 , is a function of the current density of water splitting.
For electrode with nanovoids in acidic media, a significant deviation from water equilibrium is discovered close to the electrode surface, ca. ≤ 2 μm, indicating pH is not viable to assess its impact on such electrochemical reactions involving hydroxide ions in the limited region. The concentration of hydroxide is rather close between the bottom of nanovoid and the field of the electrode, which defies the claim that pH is responsible for the damascene cobalt electrochemical plating in the nanovoid.