Electrochemical Hydrogen Evolution: H⁺ or H₂O Reduction? A Rotating Disk Electrode Study

It is a very plausible assumption that the recombination of H⁺ and OH⁻ ions in Reaction 3 is an extremely fast process and dictates instantaneous equilibrium. This fast pre-equilibrium assumption enables us to get a numerical solution of Equations 7a–7b by the application of a constraint to the concentration of H⁺ and OH⁻ ions, namely that

Here K_w = 10^{− 14} is the (by definition, dimensionless) autoprotolysis constant of water and 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 and ) 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 and 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).

• Realization of charge transfer. Charge transfer affects only the first elements of the concentration arrays; that is,
  
  
  where and .
• Realization of mass transport. The effect of mass transfer is realized in each simulation step by using the (explicit) Euler method³⁹ described by the following matrix-vector equations:
  
  and
  
  where the identity, diffusor and conveyor operators are all n × n matrices defined as:
  
  
  
• 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 and 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,⁴⁰ 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 I in both simulations, except that we assumed different values for the diffusion coefficient of hydroxide ions, . When calculating the polarization curve shown in Figure 3a and the pH profiles shown in Figure 3b we assumed that the ratio is small, having a value of 0.5, as suggested by literature.³⁰ On the other hand, the polarization curve seen in Figure 3c and the pH profiles of Figure 3d were simulated using a much larger value of β = 2 · 10⁴.

Zoom In Zoom Out Reset image size

Table I. Parameter values used in calculations at different levels of pH. Parameters marked with an asterisk are fixed; the rest of the values were determined by non-linear curve fitting (see Figure 5).

		Value at
Symbol	Definition	pH= 2.0	pH= 2.5	Unit
	bulk concentration of H+ ions	15.378	2.017	mmoldm− 3
	diffusion coefficient of H+ ions	9.171	9.804	10− 5 cm2s− 1
Eeq1	equilibrium potential, Reaction 1*	−0.341	−0.370	V
k1	reaction rate coefficient, Reaction 1	2.788	1.633	10− 7cms− 1
α1	charge transfer coefficient, Reaction 1	0.550	0.568	—
Eeq2	equilibrium potential, Reaction 2*	−1.169	−1.198	V
k2	reaction rate coefficient, Reaction 2	8.440	0.625	10− 12mols− 1cm− 2
α2	charge transfer coefficient, Reaction 2	0.598	0.861	—
ν	kinematic viscosity*	»	0.01	cm2s− 1
T	temperature*	»	298.15	K

It is apparent in Figure 3a and 3c 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³⁷:

where is defined by Equation 4 and

In the β ≈ 0.5 case, Figure 3a, a big deviation appears between the simulated polarization curve and the measured currents at electrode potentials below the limiting current region. At potentials less than − 1.4 V the reduction of water molecules already yields a significant current (Equation 5b), which at approximately − 1.6275 V surpasses the limiting current of H⁺ reduction. This potential marks a strong transition in the polarization curve: as shown in Figure 3b, at approximately − 1.6275 V the near-surface pH reaches the value of 7, the flux of H⁺ ions drops to practically zero and the total current becomes entirely dominated by the gradient of OH⁻.

The situation is somewhat different when the assumption of β ≫ 1 is applied, see Figures 3c and 3d. Now the diffusion coefficient of OH⁻ is large and this expands the effect of autoprotolysis within the solution layer. As a result, pH values at the electrode surface are nivellated (tend to be closer to neutral), however the "reaction layer"⁴¹ 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 3c 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 .

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.

The ordinary differential Equations 7a–7b 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 form

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³⁷ of the ordinary differential equation 21 is a function of the form

where and are integration constants and the Ei_s(x) exponential integral function is defined as

Then, in order to get appropriate particular solutions we formulate five boundary conditions, four fixing the values of the integration constants and a fifth defining the reaction plane position λ. The five boundary conditions are as follows:

• (i)  
  As z → ∞ the concentration of H⁺ ions tends to the bulk H⁺ concentration and the concentration of OH⁻ vanishes. Thus,

  

  and
• (ii)  
  
• (iii)  
  The gradient of the OH⁻ concentration at z → 0 is, according to Fick's first law, determined by the j₂ current density defined in Equation 5b. 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 23e, e.g., by the Newton–Raphson method³⁹ of iterative root-finding. After obtaining a numerical value of λ, the integration constants can be determined from Equations 23a–23d 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 and by the quasi-analytical approach presented here. Figure 4a 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 Figures 4b and 4c.

Zoom In Zoom Out Reset image size

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 in Equation 23e 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 , and that for the system under study (cf. to the results of digital simulation), x ≫ 1. Thus, by using the asymptotic approximation⁴⁴ that exp (x)Γ_s(x) ≈ x^{s − 1} in Equation 24, we arrive to the following approximative closed-form expression defining λ:

As shown in Figures 4b and 4c, H⁺ concentration profiles calculated from a λ value defined by Equation 26 agree well with profiles created by digital simulations assuming that β ≫ 1. Note that the analytical model yields an OH⁻ concentration of zero for the 0 ⩽ z ⩽ λ range where the OH⁻ concentration determined by digital simulation is also small. Accordingly, at the distance λ the H⁺ concentration starts to rise from a value of 0, and this is consistent with the behavior of the H⁺ profile obtained by simulation.

In Figure 4c, where the electrode potential (and thus j₂) 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 . As shown in Figure 4b, 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 while at very negative potentials j† ≈ j₂. Thus the entire polarization curve can be described as

which is a closed-form expression of the potential through Equations 4, 5b, 20 and 27, when λ is defined as in Equation 26.

The analytical expression of Equation 28 can be used quite efficiently for the estimation of parameters in non-linear curve fitting, as it is shown in Figure 5. The parameter values used for the simulation of curves are those listed in Table I. 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₁ and E^eq₂ were determined from the Nernst-equation

Here V vs. Ag|AgCl and V vs. Ag|AgCl are the standard potentials corresponding to electrode reactions 1 and 2, respectively. A constant partial pressure of bar was assumed.

Zoom In Zoom Out Reset image size

In the non-linear fitting process, the H⁺ concentration, the diffusion coefficient of H⁺ and the reaction rate parameters (k₁, α₁, k₂ and α₂) 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.³⁰ The determined kinetic parameters are in alignment with those reported earlier for H⁺ reduction on Ni electrodes.³⁵