Electrochemical Hydrogen Evolution: H + or H 2 O Reduction? A Rotating Disk Electrode Study

We study the effect of H + and OH − diffusion on the hydrogen evolution reaction in unbuffered aqueous electrolyte solutions of mildly acidic pH values. We demonstrate that the cathodic polarization curves measured on a Ni rotating disk electrode in these solutions can be modeled by assuming two irreversible reactions, the reduction of H + and that of water molecules, both following Erdey-Gr´uz–Volmer–Butler kinetics. The reduction of H + yields a transport-limited and thus, rotation rate-dependent current at not very negative potentials. At more cathodic potentials the polarization curves are dominated by the reduction of water and no mass transfer limitation seems to apply for this reaction. Although prima facie the two processes may seem to proceed independently, by the means of ﬁnite-element digital simulations we show that a strong coupling (due to the recombination of H + and OH − to water molecules) exists between them. We also develop an analytical model that can well describe polarization curves at various values of pH and rotation rates. The key indication of both models is that hydroxide ions can have an inﬁnite diffusion rate in the proximity of the electrode surface, a feature that can be explained by assuming a directed Grotthuss-like shuttling mechanism of transport.

molecules, both following Erdey-Grúz-Volmer-Butler kinetics. The reduction of H + yields a transport-limited and thus, rotation rate-dependent current at not very negative potentials. At more cathodic potentials the polarization curves are dominated by the reduction of water and no mass transfer limitation seems to apply for this reaction. Although prima facie the two processes may seem to proceed independently, by the means of finite-element digital simulations we show that a strong coupling (due to the recombination of H + and OH − to water molecules) exists between them. We also develop an analytical model that can well describe polarization curves at various values of pH and rotation rates. The key indication of both models is that hydroxide ions can have an infinite diffusion rate in the proximity of the electrode surface, a feature that can be explained by assuming a directed Grotthuss-like shuttling mechanism of transport.

Introduction
A common problem of base metal electroplating is that the deposition of electroactive metals (e.g., Zn, Co, Fe or Ni) is almost inevitably accompanied by hydrogen evolution. In aqueous solutions hydrogen evolution always occurs if current is made high enough to exceed the limiting current of the deposited metal 1 . In acidic solutions hydrogen may be formed as a result of H + reduction, while in solutions of pH > 7, the primary source of hydrogen is the electroreduction of water itself: Hydrogen evolution presents several ramifications for electrochemical plating processes. From a strictly technological point of view, the evolution of hydrogen is a serious concern because it reduces the Faradaic efficiency of metal deposition, thereby increasing the amount of time and energy utilized to deposit the desired amount of metal at a given total current density. Another problem of hydrogen codeposition is that it can affect the structural and mechanical properties of the metal and can cause embrittlement 1 .
Furthermore, hydrogen can also affect the kinetics of layer growth: if hydrogen adsorbs more efficiently on certain crystal planes of the metal, it may block these planes from further deposition so that the metal would preferentially grow on other planes 2 . This inhomogeneous growth presents difficulties for the deposition of compact metal layers at the best -that is, if the accumulated gas does not impair the entire plating process 1 .
Finally, hydrogen evolution also results in the removal of H + ions from the diffusion layer. Increasing the near-surface pH may alter the intended chemical and electrochemical reactions at the interface. For example, if the solution at the interface becomes sufficiently alkaline, it may cause the solubility of a hydroxide of the metal ion present to be exceeded, which may cause the inclusion of hydroxide or oxide species into the deposit 1,3 . This can be problematic in terms of deposit properties such as ductility and electrical resistivity. Alterna-tively, hydroxide deposition could cause the formation of relatively thick films on plated parts, thus passivating the cathode surface [3][4][5] . While the accumulation of OH − in the diffusion layer is usually negligible in sufficiently buffered solutions 6 , severe changes may occur in unbuffered or only slightly buffered From the point of view of base metal electroplating, the variation of local pH presents a challenge because it influences the kinetics of one of the most important side reactions, hydrogen evolution. However, the change of interfacial pH can, in general, affect all electrode reactions that involve H + or OH − ions in their mechanisms [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24] . Therefore, a great deal of research has aimed at the development of experimental methods for detecting local pH variations near an electrode surface. Most of these are either electrochemical or optical methods.
Optical methods usually involve the addition of a suitable pH indicator and detection by means of confocal laser scanning microscopy [13][14][15] , an approach resulting in pH values averaged over a region in front of the electrode with a thickness depending on the spatial resolution of the instrument.
A straightforward way of electrochemical pH detection is to place a pH sensitive indicator electrode close to the working electrode surface [16][17][18] : as the indicator electrode is carefully positioned in the vicinity of the cathode, its potential is changed in response to the pH of its environment. This method is capable to map the variation of pH as a function of distance to the electrode surface if a precise positioning system is used such as in scanning electrochemical microscopy 19 . Although the indicator used in this approach is usually a microelectrode, special care must be taken so that the indicator does not shield the flux of reacting species to the working electrode. In addition, the effect of IR drop must be taken into account 25,26 : that is, the method is only applicable when low cathode currents are flowing in a solution of low electrical resistance.
The problem of shielding, albeit not the problem of IR drop, can be efficiently circumvented by using the collector electrode of "generator-collector" systems in a potential sensing mode. Such experiments were described both for the rotating ring-disk electrode 20,21 and for the impinging jet cell geometries 22 . These methods are by nature not capable to "map" the pH profile near the cathode. However, by solving the diffusion-convection equations in the system, one can determine the pH in the vicinity of the working electrode based on the measured "collector" electrode potential. The detection of local pH changes immediately adjacent to the electrode surface during a redox reaction is, in addition, also possible by attaching a molecular pH probe to the cathode itself and carrying out generator-collector experiments using a single electrode 23,24 .
As can be seen above, there is a variety of experimental techniques available for determining pH changes adjacent to an electrode surface. However, there are always some experimental difficulties present and most of these techniques cannot be used under operando conditions in a metal plating cell. Nonetheless, as hydrogen evolution is an ubiquitous side-reaction of metal deposition, a selective control of its rate as well as an understanding of its effects on the technological process are vital 1,3 . Hence, there is a strong need for understanding the kinetics of hydrogen evolution in solutions of different pH.
In this paper we use numerical simulations for modelling the hydrogen evo-  (2): Taking Reaction (3) into consideration is necessary not only to determine steady-state pH profiles corresponding to given electrode potentials 29 , but also to account for the slight variation of the "limiting" H + -reduction current with the electrode potential. This latter effect often causes inaccuracies in the deter-  was applied for 30 s in order to remove any accumulated H 2 from the surface.
The current measurement was then repeated with other potential and rotation rate settings. The measured data were IR-corrected post-experimentally -the solution resistance can be determined by means of high-frequency impedance measurements.
The measurements were automated by using an Autolab PGSTAT128N potentiostat in connection with a PINE AFMSRCE rotator and by the application of the Nova v11.0 software.

8
Stationary polarization curves recorded on a Ni RDE show a strong pH dependence ( Figure 1). In case of pH 1.5 the measured current is almost entirely dominated by the reduction of H + , Reaction (1), and due to the high H + concentration there is hardly any mass transfer limitation to this electrode reaction.
On the other hand, in case of pH 3.0 (very low H + concentration) practically no H + reduction can be measured: significant currents are only yielded by the reduction of water, Reaction (2), at very negative potentials.
In case of mildly acidic pH values (2 pH 3) the afore-mentioned reactions both play a significant role in determining the shape of stationary polarization curves. From positive to negative electrode potentials, the polarization curves recorded at mildly acidic pH values start with a more-or-less exponential section where the current, due to the reduction of H + , is limited by charge transfer.
At more cathodic overpotentials the rate of charge transfer increases and the measured current is limited more and more by the transport of H + from the bulk of the electrolyte solution to the electrode interface. As mass transfer becomes the only limiting step, the stationary polarization curves exhibit an almost but not quite horizontal limiting current section. The limiting current is in reasonable agreement with predictions of the Levich equation 33 , see Figure   2: where D H + is the diffusion coefficient and c ∞ H + is the bulk concentration of H + ions, ω = 2πf is the angular frequency of rotation and ν is the kinematic viscosity of the solution. As the overpotential exceeds the limiting current section in the negative direction, the cathodic current increases again due to the reduction of water molecules, Reaction (2).
With respect to Figures 1 and 2, the following comments are of relevance: i.) The shape of the polarization curves, especially at more cathodic potentials, remain essentially unchanged when the measurement is carried out in a solution saturated with H 2 gas. This leads us to assume irreversible kinetics (i.e., the rate of Reactions (1) and (2) in the backward direction is negligible). ii.) The limiting current sections of the polarization curves recorded in solutions of mildly acidic pH values are not perfectly constant but show a slight dependence on the electrode potential (this is best shown in Figure 2).
In order to describe the measured polarization curves we shall, in what follows, attempt to develop two models that describe the kinetics of the studied electrode process where both the reduction of H + ions and the reduction of water molecules play a significant role in determining the current. As a first step we devise a numerical method for the simulation of the electrode process; then, based on the results of numerical simulations, we develop an analytical model that can explain the behaviour of the system under study and that can even be used for the fitting of measured polarization curves and for the determination of kinetic parameters.
We start with the simplifying assumption that Reactions (1) and (2) both follow Erdey-Grúz-Volmer-Butler kinetics 34 and the currents yielded by these reactions can be described as and where k i is the rate coefficient, α i is the charge transfer coefficient and E eq i is the equilibrium potential of the i th electrode reaction (Reaction (1) or (2)). In Equation (5.a), c 0 H + denotes the near-surface concentration of H + ions. Here we note that with this assumption we only wish to state that there exists an exponential dependence between the stationary partial currents and the electrode potential 35 ; in what follows we shall not draw any conclusion whatsoever about the exact mechanism of hydrogen evolution (and whether it follows a "pure" Volmer, a Volmer-Heyrovsky or a Volmer-Tafel mechanism 35 ).
In order to determine steady state currents one has to assume that the concentration of H + and OH − ions in the system is stationary, thus and Here v denotes the vector field of (stationary) fluid flow, and the symbols k +3 and k −3 were introduced to denote the forward and backward rate coefficients of Reaction (3), respectively. It was assumed at this point that mass transfer occurs only by means of diffusion and convection, while other means of transport (e.g., migration) are ignored. It was further assumed that the diffusion coefficients D H + and D OH − are constants, independent of the concentrations and of spatial coordinates.
Assuming that the radius of the disk electrode is appropriately big, the effect of transport in the radial direction can be neglected 36 and thus instead of solving the partial differential Equations (6.a)-(6.b), a solution of the set of ordinary differential equations is sufficient. Here we invoke the usual assumption 37 that in the vicinity of the electrode plane the magnitude of fluid flow in the axial direction can be approximated as where B = 0.51. Higher terms in Equation (8) Here K w = 10 −14 is the (by definition, dimensionless) autoprotolysis constant of water and c − − = 1 mol dm −3 is the standard concentration. Note that with Equation (9) we assume unit activities for the H + and OH − ions.
In the finite element simulations we tile the space under the RDE in the z direction to an n number of small volumes (cylinders of height ∆z). Accordingly, the concentrations of H + and OH − ions will be ordered into one-dimensional arrays (of size n, denoted by c H + and c OH − ) so that the i th entry of each array corresponds to the respective concentration in the i th small cylinder counted from the electrode surface.
The simulation starts from a state when both H + and OH − are evenly distributed, thus for every i and In every step of the simulation (corresponding to a time interval of ∆t) new c H + and c OH − arrays are calculated. To do this, we "update" the array entries as described in what follows (the operator ":=" in the following equations means that "updated" arrays are calculated using "present" concentration values).
1. Realization of charge transfer. Charge transfer affects only the first elements of the concentration arrays; that is, and where the I identity, D diffusor and V conveyor operators are all n × n matrices defined as: 3. Recombination of ions. In order to account for the recombination of H + and OH − ions according to Reaction (3), the following algebraic equation is solved for every i entry of the concentration arrays: Then, using the physically relevant x i roots only, and After setting the electrode potential E to the desired value, the above iteration is run until convergence is reached; i.e., until the concentration profiles become stationary. Note that step 3 assures that the product of c H + i and c OH − i entries is in accordance with the autoprotolysis constant at every i index: thus, by means of this method, one can determine a single pH profile as opposed to two "independent" profiles for H + and OH − . Current, as usual in digital simulation 40 , is calculated by approximating the near-surface concentration gradients with finite differences. This yields two partial current densities, and the sum of which gives the total current density j.
Results of two different simulations (polarization curves and pH profiles corresponding to various electrode potential values) are presented in Figure 3 for the case of pH = 2.0 and f = 900 min −1 . We used the same parameter set listed in Table 1 in both simulations, except that we assumed different values for the diffusion coefficient of hydroxide ions, D OH − . When calculating the polarization curve shown in Figure 3(a) and the pH profiles shown in Figure 3 It is apparent in Figure 3(a) and (c) that regardless of the value of β there is a good agreement between measured and simulated currents up until the point where the electrode potential exceeds the limiting current region. Currents measured at this potential range seem to be in accordance with the Koutecký-Levich-equation 37 : where j , H + is defined by Equation (4) and In the β ≈ 0.5 case, Figure 3 values at the electrode surface are nivellated (tend to be closer to neutral), however the "reaction layer" 41 is essentially broader than it was in the case of β < 1.
In case when β 1, simulated current densities are in a closer agreement with measured data and the transition of partial currents is much smoother than it was in the previous case.
It is worth to note in Figure 3(c) that the partial current calculated from the gradient of OH − already gives a significant contribution to the overall current at potentials when the rate of water reduction is still very small. This is due to the autoprotolysis process which in case of β 1 can assure that any near-surface decrease of the H + concentration would ultimately lead to an increase of c OH − . Although the numerical simulation described above is in agreement with measured data, the numerical model by nature cannot give an analytical formula that could be used for the fitting of measured polarization curves and for the quick determination of reaction parameters. Therefore, in what follows we attempt to describe the problem using an analytical approximation.

Approximate Analytical Solution of the Transport
Problem with the Constraint of Chemical Coupling The ordinary differential Equations (7.a)-(7.b) could not be solved analytically due to a strong non-linear coupling introduced by the autoprotolysis process, Equation (3). In the previous section we discussed that the upper segments of the measured polarization curves (at potentials more positive than the limiting current section) can well be described by the simple Koutecký-Levich equation (19), and we saw that problems start to arise when the flux of OH − ions also gives a significant contribution to the measured current.
In what follows we describe an essentially analytical model of the problem that tries to circumvent the non-linearity introduced by Reaction (3). We will assume that at very cathodic potentials where both H + and OH − ions have a role in determining the current density, the solution layer close to the interface can be divided into two parts. In one part, closer to the surface, only OH − ions can exist in the solution, while in the other part further away, only H + ions are present. The two regimes are separated by a "reaction plane" at a distance λ measured from the surface, acting as a sink for both ions.
This scenario, which is somewhat similar to that described by Chapman 42,43 for other reacting systems, assumes that within the two distinct regimes the concentrations of the ions obey the general convection/diffusion equation of the Here the index i stands for the species H + and OH − , and the differential equation has to be solved for both ions.
First we note that the general solution 37 of the ordinary differential equation (21) is a function of the form where χ 1,i = lim ii.) iii.) The gradient of the OH − concentration at z → 0 is, according to Fick's first law, determined by the j 2 current density defined in Equation (5.b). Thus, iv.) At z = λ, there is an equal concentration of H + and OH − ions. From this follows that v.) Finally, at z = λ, the sum of the fluxes of the two species must be zero.
This gives the following transcendental equation defining λ: Now at this point it is already possible to determine the value of λ from Equation (23.e), e.g., by the Newton-Raphson method 39 of iterative root-finding.
After obtaining a numerical value of λ, the integration constants can be determined from Equations (23.a)-(23.d) and substituted into the general solution given by Equation (22) to obtain stationary concentration profiles for both H + and OH − .
In Figure 4 we compare H + and OH − concentration profiles determined by the method of digital simulation, described in Section 4.1, and by the quasianalytical approach presented here. Figure 4(a) shows that in case of β < 1 the analytical method yields concentration profiles that do not match those obtained by digital simulation. This is due to the boundary conditions of the model which assure the equality of concentrations and the cancellation of fluxes at λ but they fail to maintain the chemical constraint dictated by Reaction (2) that at the distance λ both concentrations should have a very low value. As the results of digital simulations suggest, however, we are primarily interested in an analytical solution that is valid for the β 1 case, and this condition in itself will assure zero concentrations at the reaction plane. Thus, in a β 1 case we can expect better agreement between the simulated and analytically approximated curves, as it is also shown in Figure 4(b) and (c).
The assumption that β 1 also has implications on the analytical model, and this assumption is a step forward to obtain closed-form expressions for the concentration profiles and the current density. First, by using the substitution 23.e) and then taking the β → ∞ limit of the right-hand side, we arrive to the following expression containing λ: where Γ s (x) denotes the (upper) incomplete gamma function defined as Equation (24) is still a transcendental equation, solving it by means of an analytical approximation is however possible. First we note that Equation (24) contains the functions exp(x) and Γ 1/3 (x) with the same argument x = Bλ 3 ω 3/2 3D H − ν 1/2 , and that for the system under study (cf. to the results of digital simulation), x 1. Thus, by using the asymptotic approximation 44 that exp(x)Γ s (x) ≈ x s−1 in Equation (24), we arrive to the following approximative closed-form expression defining λ: As shown in Figure 4 In Figure 4(c), where the electrode potential (and thus j 2 ) is not very negative, the value of λ is negligibly small, and the current density, as calculated from the near-surface gradient of the H + profile, is almost exactly identical to the limiting current density j ,H + . As shown in Figure 4(b), at more negative potentials λ grows larger and the electrode is covered with a seemingly "neutral" solution layer. In this case the current density can be calculated from the gradient of the H + profile at the distance λ. This enhanced current density j † can be given as It can be shown that if the applied electrode potential is not very negative, the value of j † tends to j ,H + while at very negative potentials j † ≈ j 2 . Thus the entire polarization curve can be described as which is a closed-form expression of the potential through Equations (4), (5.b), (20) and (27), when λ is defined as in Equation (26).
The analytical expression of Equation (28) Table 1. In case of certain parameters (such as the T temperature and the ν kinematic viscosity), values were assumed and not varied during the iteration. The equilibrium potentials E eq 1 and E eq 2 were determined from the Nernst-equation Here E − −  (1) and (2), respectively. A constant partial pressure of p H2 = p − − = 1 bar was assumed.
In the non-linear fitting process, the H + concentration, the diffusion coefficient of H + and the reaction rate parameters (k 1 , α 1 , k 2 and α 2 ) were varied.
The optimized concentrations are not far from those approximated from the pH and the diffusion coefficient of H + matches its commonly accepted literature value 30 . The determined kinetic parameters are in alignment with those reported earlier for H + reduction on Ni electrodes 35 .

Discussion
In the previous section we devised two models (one based on digital simulation, another on analytical approximation) for the description of hydrogen evolution in solutions of mildly acidic pH values where both H + reduction and the reduction of water molecules play a significant role in determining the current.
Based on the analytical model it became possible to determine kinetic and transport parameters of the hydrogen evolution process by means of non-linear curve fitting ( Figure 5). The determined parameters (Table 1), when used in the simulation-based model, also yielded polarization curves that agree well with measured data (Figure 3). It is very interesting to note, however, that both models seem to "work well" (i.e., comply to measurements) only when the apparent diffusion coefficient of OH − is assumed to be very large. This property of the models deserves further explanation.
First of all we have to emphasize that D OH − in our interpretation, when used in Fick's equation of transport is indeed an apparent diffusion coefficient and thus special care must be taken when we compare it to literature values. In literature D OH − is traditionally determined by small-signal perturbations in a near-equilibrium system, e.g., by ac conductivity measurements in an alkaline solution. Implicitly utilizing the fluctuation-dissipation theorem, the aim of such experiments is to give a good measure of the diffusion coefficient as it appears, for example, in the theory of random walk. D OH − , when determined from these experiments is reasonably small (approximately one-half of the diffusion coefficient of H + ions, 45 ) but it is, again, characteristic to a near-equilibrium system.
The case studied in this paper is very different, and it is not even close to equilibrium. In our case OH − ions created at the electrode surface by the reduction of H 2 O are separated from the acidic bulk only by a very thin neutral solution layer, having a thickness not bigger than a few tens of micrometers.
Within this scheme one can assume that the driving force of transport is very big, and that the increased diffusivity of OH − can be explained by what we call a "directed Grotthuss mechanism". The Grotthuss mechanism -even in case when small-signal perturbations are applied to a system-is known to be responsible for the enhanced diffusivity of OH − compared to other ions of similar size 46 . We believe, however, that this mechanism can account for an even faster transport between the close-to-neutral electrode surface and an acidic solution region not far away from it; then, as illustrated by Figure

Conclusion
In this paper we studied hydrogen evolution as it occurs in unbuffered, mildly acidic aqueous electrolyte solutions on a Ni RDE. We demonstrated that the cathodic polarization curves measured in these solutions can be modelled by assuming two irreversible reactions, the reduction of H + and that of water molecules.
By the means of finite-element digital simulations we showed that a strong coupling (due to the recombination of H + and OH − to water molecules) exists between the two processes. We also developed an analytical model that could well describe polarization curves at various values of pH and rotation rates and that can also serve as a good basis for more detailed studies of the mechanism, e.g., by means of impedance spectroscopy.
The key indication of both models is that hydroxide ions can have an infinite diffusion rate in the proximity of the electrode surface, a feature we explained by assuming a directed Grotthuss-like shuttling mechanism of transport.   Table 1. . The two partial current densities, the catalytic (Equation (20)) and limiting (Equation (4)) current densities of H + reduction, as well as the current density yielded by the reduction of water molecules (Equation (5.b)) are also plotted, see the legend. Measured data are shown for comparison. In (b, d) stationary pH profiles are plotted for different values of the electrode potential (some are marked on the graph). The results shown in (a, b) were obtained by using a common literature value 30 of the β = D OH − /D H + ratio of 0.5 while the results (c, d) were generated by using a much higher β value of 2 · 10 4 . Other parameter values are listed in Table 1.   The applied rotation rate f is varied as shown. Dots represent measured values, curves were created by non-linear curve fitting using the Levenberg-Marquardt method 39 and Equation (28). The parameter values given in Table 1 were determined by this fitting. Figure 6: Scheme of a "directed Grotthuss mechanism" explaining the establishment of a neutral near-surface solution layer when water reduction takes place at the interface. It seems plausible to assume an extremely fast bond hopping mechanism between OH − ions, created at the interface, and the large concentration of H + ions in the acidic bulk. This mechanism can account for the large (practically infinite) apparent diffusion coefficient of OH − ions.