An analytical model of hydrogen evolution and oxidation reactions on electrodes partially covered with a catalyst

Abstract

Our previous theoretical study on the performance limits of the platinum (Pt) nanoparticle catalyst for the hydrogen evolution reaction (HER) had shown that the mass transport losses at a partially catalyst-covered planar electrode are independent of the catalyst loading. This suggests that the two-dimensional (2D) numerical model used could be simplified to a one-dimensional (1D) model to provide an easier but equally accurate description of the operation of these HER electrodes. In this article, we derive an analytical 1D model and show that it indeed gives results that are practically identical to the 2D numerical simulations. We discuss the general principles of the model and how it can be used to extend the applicability of existing electrochemical models of planar electrodes to low catalyst loadings suitable for operating photoelectrochemical devices under unconcentrated sunlight. Since the mass transport losses of the HER are often very sensitive to the H₂ concentration, we also discuss the limiting current density of the hydrogen oxidation reaction (HOR) and how it is not necessarily independent of the reaction kinetics. The results give insight into the interplay of kinetic and mass-transport limitations at HER/HOR electrodes with implications for the design of kinetic experiments and the optimization of catalyst loadings in the photoelectrochemical cells.

---

1. Introduction

This paper is a continuation of our work on theoretical simulation models used for the optimization of catalyst coatings for photoelectrochemical water splitting cells.^1,2 The model system studied is a planar TiO₂ coated silicon photoelectrode (PE) partially covered with randomly distributed Pt nanoparticles (Fig. 1). Although often criticized for scarcity and cost, Pt is still the most common catalyst used in electrochemical energy conversion devices such as in electrolyzers and fuel cells.^3–6 Employing it as finely dispersed nanoparticles helps in mitigating the material costs because it gives high electrocatalytically active surface area per gram of the material used.

Fig. 1 Illustrations of real and simulated electrodes with 50 ng cm⁻² Pt nanoparticle (d = 5 nm) loading having 2.7% surface area fraction (f_surf ≈ 0.027): (A) SEM image of Pt nanoparticles on TiO₂. Reproduced from ref. 1 with permission from The Royal Society of Chemistry; (B) a scheme of the diffusion domain approach (r ≈ 30 nm); (C) local current density distribution (mA cm⁻²) in the electrolyte for i_el ≈ −20 mA cm⁻². The grey lines are current stream lines and the black lines indicate the current density isopotential surfaces. Letters a–f indicate the boundaries in the simulation cell: a is the symmetry axis at r = 0, b the Pt surface, c the inactive electrode substrate, d the outer boundary of the cell and e indicates the geometrical position of the bulk electrolyte interface, although it is significantly farther from the electrode than the upper edge of the figure. Similarly to e, f indicates the position of the voltage contact of the electrode substrate in the 2D simulations, although in reality it is farther than the bottom of the shown portion of the simulation cell.

Optimizing catalyst loadings for photoelectrochemical (PEC) cells is quite different from other electrochemical systems in one respect: their current and power densities are relatively low because they are limited by the absorbed photon flux of sunlight. For example, a fully integrated planar PEC device, which employs equal electrode areas for both the photoconversion and electrolysis reactions, can generate at most ca. 40 mA cm⁻² when operated under unconcentrated sunlight.^7–10 At such low current densities, only a small amount of Pt is needed.^1,11–13 For example, our previous work showed that a Pt loading as low as 100 ng cm⁻² is enough to run the hydrogen evolution reaction (HER) in a PEC cell at 10 mA cm⁻² current density without inflicting more than 50 mV overpotential.¹ This loading is so low that the 5 nm Pt nanoparticles used cover only 1% of the electrode surface.

Other electrochemical devices that are usually rate-limited by the transport of chemical species (reactants and products) rather than the photon flux can be designed for much higher operating current densities. In water electrolyzers and alkaline fuel cells, for example, current densities as high as 0.2–2 A cm⁻² are achieved by employing highly porous conducting electrodes that host a large number of Pt particles, typical loadings being as high as 0.05–2 mg cm⁻².^3,5,14 In this case, the catalyst surface area exceeds the projected surface area of the electrode by more than one order of magnitude.^15,16 Optimization of such porous electrodes in terms of electrochemical performance and catalyst utilization is a classical problem of electrochemical device engineering, which requires theoretical modeling of the coupled reaction kinetics and mass transport within the porous electrode structure.^17–19 The corresponding geometrical modeling of catalyst performance in PEC cells with a low surface coverage has received less attention, perhaps because it was previously thought that no catalyst would perform well enough to be used in such low amounts.

Our first paper on this topic presented a 2D numerical model that combines the effects of mass transport and HER kinetics to describe the current–overpotential characteristics of such PEs.¹ Our second paper extended the model by considering the transport of both gaseous and dissolved H₂, as well as mass exchange between them in the electrolyte.² In this work, we present an analytical model that captures the phenomena described by the 2D numerical model and is therefore readily applicable for interpreting experimental data and performing optimization calculations in practical research. We also discuss the assumptions made with the mass transport geometry and the reaction kinetics. Although our primary interest is in the PEs, our model focuses on the mass transport in the electrolyte and the behavior of the metallic catalyst on the electrode surface. When the current flows solely through the catalyst particles, the PE–catalyst and the catalyst–electrolyte segments can be treated separately because the potential level in the catalyst can be used as a boundary condition for both segments. To sum up, the potential differences that correspond to the same current density and catalyst potential yield the total voltage, as with fully covered PEs.²⁰ This way, the model developed in this paper can be connected to models describing the photovoltaic operation of the PE in a complete PEC device.

The paper is organized as follows. We first present the main observations from the previous two papers with regard to the geometrical considerations to give a rationale for the 1D modeling approach in this paper (Section 2). Thereafter, we discuss our general approach and the details of the kinetic model coupled with the mass transport, and how we solve the equations. As the mass transport is included in our model, the effects of the limiting current density on reaction kinetics are discussed, as well as how the reaction mechanism affects the limiting current density. The relative simplicity of our model allows us to derive analytical expressions that describe the micropolarization region and high overpotentials (i.e. “Tafel equation”), which are compared with the full numerical solution.

2. Current distribution and mass transport losses at the nanoparticle electrode with a low surface coverage

One of the observations from earlier 2D numerical simulations of our model system (Fig. 1) was that when the current density per electrode area was considered, the mass transport losses were almost independent of the Pt loading,^1,2 which allowed simplifying the mass transport problem to one dimension. As a result, an analytical 1D solution, mathematically similar to the diffusion of only dissolved H₂ but including the effects of the bubble transport and dissolution kinetics, was found and shown to be valid at current densities below ca. 40 mA cm⁻².² Why the catalyst loading had only a little effect on the mass transport could be rationalized by realizing that the stagnant diffusion layer was much thicker than the average distance between the Pt particles. Because the flux profile is a result of a linear drift–diffusion model, it can also be considered as a superposition of the overlapping, hemispherical flux profiles of the individual catalyst particles:^21–23 Close to the electrode surface the flux is mainly affected by the nearest particle, creating the hemispherical profile, whereas far from the interface the effects of the individual particles are blurred, forming a uniform 1D profile.²³ Numerical simulations (Fig. 1C) illustrate this by showing how the 1D flux converges to a spherical flux towards the catalyst particle only very close to the electrode surface. The fact that going from a uniform Pt coating to sparsely distributed Pt nanoparticles hardly affects the mass transport losses indicates that the convergence of the transport flux to the Pt particle brings negligible contribution to the total mass transport resistance. The 1D mass transport picture holds as long as the average distance between the nearest particles is much shorter than the diffusion layer thickness that corresponds to the limiting current density of the mass transport.^23,24 Therefore, the 1D picture does not apply to the extremely sparse catalyst decorated electrodes and to very rapid mass transport. Since at least the sparse catalyst array could be a relevant scenario, in addition to other study topics, we also aim to identify when the 1D picture is no longer accurate.

Another key observation was that the mass transport and kinetic overpotentials depend on different current densities: the mass transport overpotential depends on the current per electrode surface area whereas the kinetic overpotential arises from the current density per catalyst surface area, which in our case is smaller than the electrode area. This separation suggests that our 2D model could be fully simplified to 1D by coupling these current density scales to each other. In this work, we show that this is indeed correct: the 1D analytical model obtained agrees remarkably well with the 2D numerical simulations. This result suggests that other 1D models or point models of the homogeneous electrode–electrolyte interface could also be extended to electrodes partially covered with a catalyst.

3. One-dimensional model of HER/HOR on catalyst particle arrays

3.1. General modeling principle

We describe the electrode using a 1D model. The concentrations on the catalyst surface are determined by the mass transport model and the current density per electrode area (i_el), whereas the total overpotential is given by the model of the reaction kinetics, current density per catalyst area (i_cat) and surface concentrations. These current density scales are coupled together via the ratio of the catalyst surface area (A_cat) to the geometrical electrode area (A_el), f_surf

i_el = f_surfi_cat (1a)

In the context of the present paper, eqn (1) is the general principle for coupling a kinetic model that considers varying catalyst surface areas to a mass transport model that neglects them.

In a broader context we point out that, mathematically speaking, f_surf can describe equally well also by varying the catalytic activity of the electrode surface because in all the equations of our model, it always appears as a factor multiplying the exchange current density of the reaction. By scaling the exchange current density of the reaction with the surface area fraction, we assume that the catalyst distribution on the surface is sufficiently homogeneous such that the changes in the catalyst loading can be described as changes in the catalytic activity of the electrode, without considering the mass transport in the vicinity of the particles.

The conventional way of using f_surf is to describe electrodes with an active area that is larger than the geometrical electrode area (f_surf > 1, Fig. 2C), i.e. as the surface roughness factor.^15,25,26 Using the ratio to describe surfaces that are only partially covered with the catalyst (f_surf < 1, Fig. 2A) is apparently not as common, and we have found only a few examples of this.^24,27 Our kinetic model is based on earlier models,^15,28 and the details of the reaction kinetics are discussed in several earlier articles.^15,28–34 In this article, we verify that eqn (1) is as accurate as our 2D model of nanoparticle arrays and thereafter analyze the effects of the catalyst surface area on the electrode operation. Although our discussion does not go into the details of the reaction kinetics, such as the effects of the kinetic parameters^28,34,35 or blocked adsorption sites,¹⁵ the results from previous theoretical analyses can be applied to our situation, namely the study of the effect of the catalyst surface area by varying either the exchange current density or the mass transport limitation.

Fig. 2 Cross-sectional schemes of the different regimes of the f_surf values: (A) partially covered electrode, (B) smooth, planar electrode and (C) rough, planar electrode. The grey color indicates the electrochemically active material.

3.2. Mass transport and surface concentrations

We assume linear mass transport to allow direct comparison with our earlier 2D simulations.¹ For the sake of simplicity and to allow analytical solutions to be reached, the effect of H₂ bubbles on the transport of the dissolved H₂ molecules is neglected here.²

In practice, linear mass transport means that the surface concentrations can be described using a bulk concentration and a limiting current density (which depends on the bulk concentration)

The concentration of a species n on the catalyst surface is denoted by c_n, the bulk concentration is indicated by superscript b and the concentration that corresponds to the thermodynamic reference potential (in this case the standard hydrogen electrode, SHE and V_SHE) by superscript 0. The SHE conditions correspond to 1 mol l⁻¹ (1 M) proton concentration and 1 bar H₂ pressure, i.e. 0.77 × 10⁻³ M, but it is possible that only a small fraction of H₂ is in a form that can participate in the reaction at the electrode surface.³⁶ In our simulations, c^b = c⁰ and the bulk concentration is included for the sake of generality.

For SHE conditions, we use the limiting current densities that we determined from earlier 1D simulations (“flat Pt” in ref. 1, for 5 μm diffusion layer i_lim,H⁺ = −35 940 mA cm⁻² and i_lim,H2 = 15.4 mA cm⁻², see Table 1 for the simulation parameters). We follow the convention that positive currents correspond to HOR and negative currents to HER.

Table 1 The simulation parameters and their original source (if not based only on 1)

Symbol	Explanation	Value
T	Temperature	298.15 K
i lim,H2	Limiting current density of H2 transport	15.4 mA cm^−2
i lim,H^+	Limiting current density of H^+ transport	−35 940 mA cm^−2
	Equilibrium proton concentration	1.0 M
	Equilibrium H2 concentration	0.7698 mM^36
i 0,V	Volmer reaction exchange current density	100 mA cm^−2^16
r H	Ratio of Heyrovsky and Volmer rates	0^16
r T	Ratio of Tafel and Volmer rates	9.5
θ ^0	Equilibrium H-coverage	0.67^46

---

In addition to simplifying the details of the catalyst array geometry and mass transport to the single parameter f_surf (Sections 2 and 3.1), we focus on the steady state operation, and therefore omit the transient behavior. The phenomena related to the electrode geometry that f_surf alone cannot describe, such as the surface concentration transient in a nanoelectrode array,^22,23 are out of the scope of our model. Although surface diffusion and the spillover effect may be important with catalyst arrays, especially with rapid mass transport,²⁷ we have neglected them here, as it seems that the HER/HOR on Pt could be described only by mass transport and reaction kinetics.²⁴

3.3. Reaction kinetics

The total hydrogen oxidation and evolution reaction proceeds through intermediate reactions called Volmer, Tafel and Heyrovsky steps

H₂ ⇄ 2H⁺ + 2e⁻ HER/HOR (3a)

H₂ ⇄ 2H_ads Tafel (3b)

H₂ ⇄ H_ads + H⁺ + e⁻ Heyrovsky (3c)

H_ads ⇄ H⁺ + e⁻ Volmer (3d)

Hydrogen atoms adsorbed on the surface are marked by H_ads and the reaction proceeding from left to right corresponds to the HOR (and from right to the left to the HER). Although the HER/HOR on Pt has been studied extensively, even the most recent literature is somewhat inconclusive about the reaction mechanism, with some results indicating the Volmer–Tafel (V–T),^16,37,38 and others the Volmer–Heyrovsky (V–H) mechanism³⁹ as the dominant path. To allow direct comparison with our earlier 2D simulations¹ where we assumed the V–T path with Volmer as the rate determining step (RDS), we assume the same mechanism in this article as well. Some features in our analysis are specific to the V–T mechanism, but the general principle of the model and the applicability of the 1D model should be independent of the reaction mechanism.

The V–T mechanism is convenient for the analysis of the model and simulations, because we can derive an analytical expression for the fraction of the reaction sites covered with H_ads(θ) as a function of the current density. Alternatively, we can also determine the current density as a function of θ, as discussed in the next section. We begin with the reaction rates of the Tafel (ν_T) and Volmer (ν_V) steps

The rate of the Volmer step (eqn (4b)) is described with an extended version of the Butler–Volmer (B–V) equation that takes into account the proton concentration and the fraction of reaction sites occupied by adsorbed hydrogen atoms. The elementary charge is denoted by q_e, the Boltzmann constant by k_B and temperature by T (in Kelvin). The symmetry factor of the Volmer step is assumed to be β = 1/2 (and thus η′/2 in the exponent), as is common. The superscript 0 marks the equilibrium value at the thermodynamic reference potential (i.e. the SHE potential) also in the case of θ. The exchange rates of the reactions are denoted by subscript 0 and the equilibrium potential of the reaction is V₀ (0 V vs. SHE) with V being the potential of the catalyst surface. Because the resistive losses in metallic catalysts are negligible, the potential of the catalyst particles can be considered to be constant, and therefore, the catalyst potential V can be used directly as a boundary condition to connect a photovoltaic model of the PE to the kinetic model of the catalyst. This works as long as the HER/HOR proceeds only on the catalyst surface and hence all current flows through the PE–catalyst contact and none, or a negligibly small fraction, through the PE–electrolyte interface. Note that this also holds in the case of pinched-off catalyst particles⁴⁰ when the charge transport to the PE–catalyst contact depends on the nearby PE–electrolyte interface, because, even in this case, the catalyst potential can be used as the reference point for the calculations. On the other hand, if the reaction kinetics depends not only on the catalyst particle and its potential, our approach of separating the catalyst operation from the rest of the electrode does not apply, and another modeling approach should be adopted instead. This would be true for instance if the catalytic activity of the bare PE surface could not be neglected, or the metal coverage was thin enough to allow the substrate to affect the electronic properties of the catalyst surface, making it different from the surface of the bulk material, thus affecting the reaction kinetics.^41,42

Eqn (4) gives the reaction rates on the catalyst surface. We substituted current density for the surface concentrations in eqn (4a) and (4b) to explicitly show the variables that we use in the calculations. Because no electrons cross the electrode–electrolyte interface in the Tafel reaction, the Volmer rate multiplied by the Faraday constant and f_surf yields the current density per electrode area.

The exchange current density is i_0,V = F·ν_0,V, where F is the Faraday constant. When the current density and θ are known (see Section 3.5), the overpotential is determined numerically from this equation. In addition, it is possible to derive analytical expressions for the micropolarization range and “Tafel equations” for high negative and positive overpotentials (Section 3.7).

The overpotential in eqn (4c) is the total overpotential and includes the losses from both the reaction kinetics and mass transport. The mass transport overpotential is the Nernst potential difference between the catalyst surface and the bulk electrolyte.

The current density dependence can be determined by taking logarithms of eqn (2a) and (2b). Often, it can be assumed that the limiting current density of the proton transport is significantly higher than the current density (|i_el/i_lim,H⁺| → 0) or the H₂ mass transport limitation, leading to a solution that depends only on the limiting current density of the H₂ transport^1,16,43,44 This reflects the starting point of our model. When the current density per electrode area is considered, the mass transport losses are independent of the amount of catalyst (f_surf). If the H₂ mass transport limitation is lower than the exchange current density of the electrode, then not only is the HOR under diffusion control, but also the HER overpotential is approximately equal to the mass transport overpotential, thus effectively determined by the H₂ transport.^16,43

In equilibrium, the net rate of both reaction steps must be zero, from which we can derive the potential for the reversible hydrogen electrode (RHE) with the given bulk concentrations. Starting from eqn (4), this condition yields an expression similar to eqn (6) (but with the logarithm of instead of ), so the kinetics of our model is consistent with the thermodynamics of the reaction. It would also be possible to derive expressions for and numerically determine the dependence of the exchange current density with respect to the proton and H₂ concentrations, but this property and comparison to experimental results^37–39 are beyond the scope of our analysis.

3.4. Simulation parameters

Table 1 contains the simulation parameters introduced in the previous sections. We used the same values that we used previously¹ or the values derived from these simulations. The limiting current densities correspond to the solution of the 1D drift–diffusion problem over a 5 μm thick diffusion layer. In the case of protons, electroneutrality is enforced and perchlorate ions (ClO₄⁻) are the anions in the electrolyte.

In the 2D simulations, we assumed spherical particles that were simulated in cylindrical geometry. Therefore, the catalyst surface area fraction is

where r is the radius of the catalyst particles, f_exp the fraction of the catalyst particle surface that is exposed to the electrolyte and R_cell the radius of the simulation cell. Previously, we assumed spherical particles (5 nm diameter, i.e. r = 2.5 nm) that were slightly embedded in the substrate, so that 95% of their surface area was exposed to the electrolyte (i.e. f_exp = 0.95).¹

The radius of the simulation cell corresponded to a given catalyst mass loading (ng cm⁻², L_cat)

The density of the catalyst material is ρ_cat, which for Pt is ρ = 21 450 kg m⁻³.⁴⁵ Depending on the regularity of the particle placement on real electrodes, eqn (9) may slightly overestimate either the center to center distance between particles (d_c–c = 2R_cell) if starting from the mass loading or the catalyst loading if starting from the average distance between the neighboring particles. In both cases, the overestimation originates from the use of cylindrical/circular unit cells. In regular particle arrays, the unit cells are regular polygons and the distance between the particles is twice the apothem (the distance between the center of the polygon and the midpoint of its side). Thus, for a given unit cell area (particle density) the radius is larger than the apothem and for a given interparticle distance, the area of the circle is smaller than the area of the polygon.

3.5. Current density and hydrogen coverage of the catalyst

As mentioned, eqn (5) includes three unknown variables, i_el, θ and η′. We can use either the current density or the hydrogen coverage as the starting point and calculate the other from the Tafel rate. In steady state, θ does not change over time, yielding The factor 2 comes from the stoichiometry of the reactions (two hydrogen atoms per Tafel step, but one per Volmer step). We also define r_T = ν_0,T/ν_0,V (similarly to Wang et al.¹⁵), the relative exchange rate of the Tafel step compared with the Volmer step. Although one can consider r_T > 1 to correspond to the Volmer step as the RDS and r_T < 1 to the Tafel step as the RDS, the value of r_T needs to be significantly lower than 1 to produce the 30 mV per decade Tafel slope that is considered an indication of the Tafel step being the RDS.³⁰

Since the only unknown variables in eqn (10) are the current density and the surface hydrogen coverage, we can use either of them as the starting point and determine the other analytically. In our case of the V–T mechanism and possibly also in the case of other Langmuir–Hinshelwood reactions, θ is a mathematically convenient starting point. However, because the convenience may be limited to the specific combination of the kinetic and mass transport models and because θ is likely a less intuitive starting point than the current density, we also discuss our model with the current density as the starting point.

Calculating the current density (i.e. starting from θ ∈ [0,1]) obtained using eqn (10) gives

The numerator is the Tafel rate with the bulk H₂ concentration and the denominator corresponds to the effects of the H₂ transport. It is readily apparent that the proton concentration has no effect on the interdependence of the current density and θ. Using eqn (11) it can be shown that ∂i_el/∂θ < 0 (when θ⁰ ∈ (0,1) and θ ∈ (0,1)), so the extremes θ = 0 and θ = 1 correspond to the kinetic limits of the HOR and HER current densities, respectively. Note that in the case of the HOR the mass transport limitation remains and affects the limiting current density. This is discussed in more detail in Section 3.6.

To determine the hydrogen coverage from the current density, eqn (10) is reorganized to the standard form

Aθ² + Bθ + C = 0 (12a)

Although there are minor differences, this is equivalent to the solution of Wang and co-workers.^15,47 The physically sensible solution corresponds to the negative root, because it yields θ = θ⁰, when i_el = 0 (and ). The last term of C corresponds to the kinetic limiting current density: when we assume that either the surface is fully covered with hydrogen or the coverage tends to zero (and ), the highest current densities per catalyst area that the Tafel step allows are Therefore, the last term can be rewritten as This term is therefore the mathematical way to express the effect of the finiteness of the Tafel exchange rate on the hydrogen coverage as a function of the current density. We denote the kinetic current density limitations per electrode area in a shorter form that includes the effect of the catalyst surface area

i_lim,HER,kin = f_surfi_{lim,HER,kin,Pt} (14a)

i_lim,HOR,kin = f_surfi_{lim,HOR,kin,Pt} (14b)

In the micropolarization region, the current density is small compared with the exchange current density of the reaction. Regardless of which step is the RDS, the last term of the C is at most equal to the current density divided by the effective exchange current density of the electrode. If the Volmer step is the RDS, r_T > 1 and the term is even smaller than this fraction. On the other hand, if r_T < 1, the denominator is equal to the exchange current density of the reaction, so the assumption is valid irrespective of the value of r_T. If the last term of C can be neglected, eqn (12c) is simplified to Note that at low current densities, i.e. i_el/i_lim,H2 → 0, this expression becomes constant and depends only on the equilibrium coverage and the bulk H₂ concentration, which implies that if the bulk H₂ concentration is fixed by the experimental conditions, the H₂ transport determines the current density range, where the constant hydrogen coverage can be assumed in the interpretation of the experimental results. Although the hydrogen coverage is independent of f_surf, low catalyst loadings allow the current density per catalyst area to reach higher values before the hydrogen coverage changes significantly from its equilibrium value.

3.6. Current limitation due to kinetics, mass transport and their combination

As discussed briefly in the previous section, in addition to the mass transport, the reaction kinetics also limits the current density in the case of the V–T mechanism because the rate of the Tafel step depends directly on θ, but not on the overpotential. Therefore, the steady state current density cannot be enhanced beyond the limits imposed by the Tafel step. The expressions for this limit for both HER and HOR (under the SHE conditions) are given in eqn (13).

Because proton transport does not affect θ directly, we can consider the mass transport and kinetic limitations of HER separately. Therefore, if the limiting current density of proton transport is lower than the limit of the reaction kinetics, the catalyst is never fully covered with adsorbed hydrogen atoms, and the current is limited by the proton transport. In the opposite case, the current density is limited by the Tafel step (eqn (11) and (12), θ → 1 corresponds to i_el → i_lim,HER,kin) and the mass transport limit is not reached. In eqn (15) the kinetic limitation is neglected and therefore θ → 1 corresponds to i_el → −∞.

Because the mass transport losses are mainly dictated by the H₂ concentration (eqn (6) and (7)), the limiting current density of the H₂ transport is an important parameter, even when considering only the HER.^1,16,44 However, as already mentioned, the reaction kinetics also limits the HOR current density. Therefore, although it may be easy to measure the HOR limiting current density, it might not be equal to the H₂ transport limitation. Moreover, because the H₂ transport and θ are directly coupled to each other in the Tafel rate, mass transport and kinetic limitations cannot always be considered separately, like in the case of the HER. As mentioned with regard to eqn (11), the current density is maximized at the limit θ → 0, so the corresponding current density is the HOR limiting current density.

We recognize the kinetics-limited HOR current density from eqn (13b) and (14b) and (after multiplying both the numerator and denominator by i_lim,H2) rewrite this as Setting θ = 0 in eqn (12a), i.e. C = 0, yields the same result. (The other possibility that corresponds to C = 0, θ = −B/A always yields a coverage that is either negative or higher than 1, and therefore unphysical.) The terms inside the parentheses in eqn (16b) correspond to the two reasons for the catalyst surface to be devoid of hydrogen: at currents near the mass transport limit (first term) there is no H₂ in the electrolyte that could be adsorbed on the catalyst, whereas the kinetic limiting current density (second term, see also eqn (14b)) corresponds to the rate at which the H₂ molecules in the electrolyte are adsorbed and dissociated to hydrogen atoms on an empty catalyst surface. Note that the mass transport losses are determined by the mass transport limitation (eqn (7)), and not the total limiting current density.

From eqn (16b) we observe that, when considering the HOR limiting current density, the mass transport and kinetic limits can be considered independently only when either limiting current density is significantly higher than the other and the HOR limiting current density is (approximately) equal to the lower limitation. This dependence of the total limiting current density on the reaction kinetics and mass transport is analogous to the Koutecky–Levich equation²⁶ and it has been observed in both theoretical and experimental studies (e.g.ref. 28, 29, 32, 33 and 48). However, eqn (16b) is limited to the HOR limiting current density of the V–T mechanism and is not applicable to the limiting current densities in general. In addition, unlike in the typical application of the Koutecky–Levich analysis, we did not neglect the effect of the reverse reaction or the hydrogen coverage in the derivation of eqn (16b), so the errors associated with these simplifications do not affect our results.^30,49

The total limiting current densities calculated using the simulation parameters in Table 1 are shown in Fig. 3. At low catalyst surface ratios, the electrode is limited by the reaction kinetics and at sufficiently high loadings by mass transport. In addition, there is also an intermediate region (in our case, f_surf values range between 10⁻⁴ and 0.01) where the total current density limitation is noticeably lower than either limit alone. The largest difference occurs when the kinetic and mass transport limitations are approximately equal to each other (i_lim,HOR,kin ≈ i_lim,H2), in which case the combined limiting current density i_lim,HOR,VT is about one half of the individual limits. The difference between Fig. 3A and B is the current density normalization. When the current density per electrode area is considered (A), the mass transport limitation is independent of the amount of catalyst, whereas in the case of the current density per catalyst area (B), the kinetic limit is constant and the mass transport limit depends on f_surf (as indicated by eqn (1)). Note that because in the model the mass transport limit arises purely from the 1D mass transport, which has a constant limit per electrode surface area (Fig. 3A), the mass transport limit calculated per catalyst surface area (Fig. 3B) tends to infinity as f_surf, and thereby the catalyst area, approaches zero (i.e. constant divided by zero). However, as this happens, the average distance between the catalyst particles also becomes large and eventually similar to the diffusion layer thickness, breaking the assumption that the transport can be considered mainly as one-dimensional. This leads to a situation where the transport limitation may be considered as a combination of spherical and 1D mass transport limits, which could possibly be expressed in the same fashion as eqn (16b), as done by Montero et al.²⁴.

Fig. 3 The limiting current densities in eqn (16b) as a function of f_surf in (A) current density per electrode area and (B) current density per catalyst area. Note the different vertical axis scales. The simulation parameters are given in Table 1. The filled circles mark the points where the current limitation undergoes transitions from the kinetic control (red, f_surf = 10⁻⁴) to the mass transport control (blue, f_surf = 10⁻²), and the intermediate region where the current is limited both by the kinetics and mass transport (black, f_surf = 10⁻³). These three characteristic cases can be identified in the other figures as red, blue and black curves, respectively (Fig. 5–8).

3.7. Simplified, analytical solutions

So far, we have discussed the interdependence of the current density and the hydrogen coverage of the catalyst together with the limiting current densities. As mentioned in Section 3.3, we determine the overpotential from the Volmer current density (eqn (5)) after obtaining the current density and the hydrogen coverage. Although the exact solution must be obtained numerically, we can derive analytical expressions for the overpotential in the micropolarization region and “Tafel equations” for high negative and positive current densities, similar to the corresponding expressions derived from the B–V equation (Chapter 3.4.3 in ref. 26). For the micropolarization range this gives The Tafel equations for high negative and positive overpotentials are, respectively The last expression (eqn (17c)) is not needed, if f_surf and i_0,V are high enough compared with i_lim,HOR,VT that the limiting current density is already achieved in the micropolarization region, which eqn (17a) describes accurately (Section 4.2). As discussed in Section 3.3, θ is a function of the current density (eqn (12) and (15)) or current density is a function of θ although neither of the possibilities is expressed in eqn (17). In all our calculations, we have used either eqn (11) or (12); thus the effect of the reaction kinetics is not neglected, as in eqn (15). In the general form, without writing out the current density dependence on the surface concentrations (or any other dependencies of the physical quantities and simulation parameters), eqn (17a) and (17b) would be If we can assume that θ ≈ θ⁰ and , HER/HOR is described accurately by the concentration-independent B–V equation, and these expressions are simplified to the more familiar forms that depend only on the current density and the exchange current density (including f_surf). Eqn (15) suggests that in the micropolarization range the hydrogen coverage is approximately constant, if the current density remains small compared with the H₂ transport limitation, so the micropolarization range can probably be described using a simplified expression (the first term of the numerator in eqn (17a)), if the electrode is not mass transport limited. In addition, although in general the Tafel slope of the HER overpotential is not necessarily 120 mV per decade (at room temperature, i.e. η′/2 in the exponent in the B–V equation, as in eqn (4b)), it may be approximately equal to that up to some overpotential that depends on how rapid the mass transport is (both protons and H₂ through the surface hydrogen coverage).

4. Results and discussion

4.1. Comparison with 2D simulations

As can be seen in Fig. 4, the analytical 1D model described in Section 3 matches accurately the 2D simulations using the model introduced in ref. 1, for the whole catalyst loading range from 0.01 ng cm⁻² (f_surf ≈ 5.3 × 10⁻⁶) to 1000 ng cm⁻² (f_surf ≈ 0.53), except for the high HOR overpotentials and low catalyst loadings (1.0 ng cm⁻² or less). The match between the HER overpotentials is excellent, but at high overpotentials, the HOR limiting current densities differ, because at catalyst loadings below 1.0 ng cm⁻² (f_surf ≈ 5.3 × 10⁻⁴) the distance between the particles and its effect on the mass transport cannot be neglected and therefore the 1D approximation is no longer accurate. Although at the lowest f_surf values the 1D simulations differ from the 2D simulations due to inaccurate treatment of mass transport, the reduction in the HOR limiting current density, when f_surf is reduced, is due to the reduced kinetic limiting current density (Fig. 3A and eqn (16b)), and not mass transport.

Fig. 4 The current density–overpotential curves (iη curves) for different Pt loadings simulated using the 1D model (lines, eqn (5) and (11)) compared with the 2D simulations (markers).¹ The colors indicate the different Pt loadings.

Assuming that the mass transport behaves as described in the model in ref. 24, the spherical diffusion becomes a significant factor, when its limiting current density becomes similar to the 1D limitation. The HOR mass transport limit for a hemispherical particle can be calculated from the well-known formula.^23,26 For a particle with a 5 nm diameter (D_H2 = 5.11 × 10⁻⁵ cm² s⁻¹ and c_H2 = 0.77 mM) this limit is about 3.1 × 10⁴ mA cm⁻² per catalyst area. Because in our case almost the entire particle is exposed, the mass transport limitation per catalyst surface area would be about half of this, which coincidentally is approximately equal to the kinetic limitation with our simulation parameters (Fig. 3B). Therefore, the 1D and hemispherical mass transport limitations would be approximately equal when f_surf ≈ 10⁻³ (Fig. 3B). Considering that the HOR limiting current density of the 1D model becomes inaccurate between 1.0 ng cm⁻² (f_surf ≈ 5.3 × 10⁻⁴) and 10 ng cm⁻² (f_surf ≈ 5.3 × 10⁻³, Fig. 4), the comparison of the hemispherical mass transport limitation to the 1D limitation is probably a good criterion for estimating the lower limit of the f_surf range where the 1D model is accurate, although the size of the spherical mass transport region (and therefore its limiting current density²⁴) depends on f_surf. Because the limiting current density of the spherical mass transport depends on the particle size, the f_surf range where our model is valid also depends on it, with smaller particles allowing accurate simulations about the lower catalyst surface area fractions. Intuitively, this can be understood so that with smaller particles, a given catalyst surface area is divided among more particles, creating a more homogeneous catalyst distribution. At the other extreme, if the same catalyst surface area corresponds to a single larger particle on the electrode surface, certainly mass transport must be considered in detail.

Note that the point where the 1D model becomes inaccurate depends on the 1D mass transport limitation, lower limiting current densities (i.e. thicker diffusion layer) being accurate at lower f_surf values. The reaction kinetics also affects the accuracy of the 1D model: in general, slower reaction kinetics (lower exchange current density) corresponds to a lower kinetic limiting current density. Therefore, when the reaction kinetics is sufficiently slow, the mass transport limitation can be neglected entirely (Fig. 3B and eqn (16b)). For the HER/HOR, Pt in acid is therefore an extreme case, because the exchange current density, and thus probably also the kinetic limitation is the highest among the known catalyst materials, and therefore the results are more sensitive to the inaccuracies of the mass transport model than those with any other catalyst.

The excellent match between the 1D and 2D model means that the 1D model is accurate enough to describe the steady state behavior of these nanoparticle electrodes, especially for the HER. The match also implies that it should be possible to accurately extend 1D models derived originally for planar catalyst electrodes to describe also electrodes partially covered with catalyst particles, by scaling the exchange current density of these models using the catalyst surface area fraction that depends on the catalyst loading. Because of the excellent agreement in Fig. 4, we use only the 1D model in the rest of the paper.

4.2. The micropolarization range and the Tafel equations

Fig. 5 compares the simplified solutions for the micropolarization range and Tafel equation (eqn (17)), with the full 1D solution of the current density–overpotential curve (iη curve). In all cases, the micropolarization curve is accurate when |η| < 50 mV, and the Tafel equations are accurate for |η| > 100 mV. It is not evident from Fig. 5, but the Tafel slope of the HER branch at high overpotentials (when current is not yet limited by the kinetics or proton transport) increases as the current density increases. For example, in the case of f_surf = 0.01 the slope is approximately 120 mV per decade up to about −0.12 V vs. RHE and at about −0.15 V vs. RHE it is approximately 150 mV per decade. Since the hydrogen coverage decreases with increase in the current density, there are less free sites for proton adsorption when the HER rate increases. (The HER corresponds to the negative currents and HOR to positive currents.) Therefore (the absolute value of) the overpotential for a given current is higher than a coverage-independent B–V equation would predict, corresponding to an increased Tafel slope. Ultimately, the current density may become limited by the reaction kinetics rather than proton transport (Section 3.6), as indicated by the f_surf values 10⁻⁴–0.1.

Fig. 5 Comparison of the expression for the micropolarization range (μ, solid lines, eqn (17a)) and Tafel equations (dashed lines, eqn (17b) and (17c)) with the full 1D solution (markers, eqn (5)). The colors of the lines and markers indicate the value of f_surf.

4.3. HOR limiting current density

As eqn (7) illustrates, the H₂ transport limitation is a crucial parameter also when considering the HER, because it is the main descriptor of the mass transport losses. Since in the case of the V–T mechanism, the HOR limiting current density does not necessarily correspond to the H₂ transport limitation, it may be possible to wrongly estimate the mass transport limitation by assuming it to be equal to the measured HOR limiting current density. This could lead to erroneous estimates of the mass transport losses and kinetic overpotential. We therefore discuss this detail of the reaction kinetics and the effect of the catalyst loading on it in the following sections.

Fig. 6A shows that the hydrogen coverage decreases sharply as the current density (i_el) approaches the value of eqn (16) (i_lim,HOR,VT). This is true at both low (f_surf < 10⁻⁴) and high (f_surf > 10⁻²) catalyst surface area ratios. Most importantly, the predicted limiting current density also holds in the intermediate range of f_surf values, where the kinetic and mass transport limitations are similar to each other (i_lim,HOR,kin ≈ i_lim,H2) and significantly higher than the combined limit i_lim,HOR,VT (Fig. 3). For example, in the case f_surf = 10⁻³ the mass transport limit has the value 1.87 and the kinetic limit value 2.15 in the normalized current density scale of Fig. 6. The fact that the overpotential increases rapidly at this limit, even in the intermediate case f_surf = 10⁻³ (Fig. 6B), confirms that the HOR current density is limited by the combination of the reaction kinetics and mass transport.

Fig. 6 (A) The hydrogen coverage of the catalyst surface for different values of f_surf as a function of the current density (eqn (11)) normalized using the total HOR limiting current density (eqn (16b)). (B) The normalized iη curves for the same f_surf values are as shown in A. Note that in A, values 0.01–10 and in B, 1.0 and 10 overlap, hence they are not easily discernible.

When the current density per catalyst area (i_cat) is considered, the 1D mass transport limitation is divided by f_surf, hence the limiting current density in the kinetics-limited low f_surf cases is significantly higher than in the mass transport limited cases (Fig. 3B). Therefore, contrary to what Fig. 6A may initially seem to allude to, the hydrogen coverage decreases slowly when f_surf is decreased, as also eqn (15) suggests, so that the related discussion about the effect of f_surf is not contradicted.

Fig. 6B shows the normalized current density as a function of the overpotential. The curves have shapes similar to what would be expected from the combination of the classical B–V kinetics at low current densities and a sharp mass transport limitation at high current densities (Section 3.3 in ref. 26). This is well in line with Fig. 3 and 6A and the theoretical description behind them. As the catalyst surface fraction increases, the iη curve becomes less ‘S-shaped’, signifying the transition from kinetic to mass transport dominance (Fig. 3). As the current density (i_el) tends to the limiting current density i_lim,HOR,VT, the overpotential increases sharply due to the decrease in the surface hydrogen coverage (eqn (5)). The similarity of the shapes to the simple B–V kinetics with mass transport limitation means that it could be difficult to distinguish from experimental curves, whether their current limitation is due to mass transport or kinetics, unless the catalyst loading (surface ratio) was varied systematically in the experiment. The mass transport limited case can probably be recognized, because its iη curve is determined by the mass transport limit alone,^1,16,43,44 but other limiting current densities are likely more difficult to determine without data about several catalyst loadings.

4.4. Validity of the model in the general Volmer–Heyrovsky–Tafel case

So far, we have considered only the V–T path, neglecting the Heyrovsky step. Although near equilibrium the current density of the reaction dominated by the V–T path could be described accurately without considering the V–H path, it likely affects the current density, especially at high overpotentials.^15,32,47 Unlike the Tafel rate, the Heyrovsky rate is directly affected by the applied potential; thus its rate can be increased with the overpotential.^15,28 Therefore, the total current density is limited only by the mass transport limitation, allowing the current density to exceed the V–T limit, which would be observed as an additional intermediate shoulder or plateau in the current iη curve.^15,28

Such shoulders in the HOR current density curve have been observed experimentally at least by Chen and Kucernak⁴⁸ and Elbert et al.;⁴⁷ Wang et al.¹⁵ interpreted the experimental results of Chen and Kucernak⁴⁸ in this way. Although our planar electrode model differs geometrically from the single Pt particle measurement setup,^48,50 the effects of the particle size in the measurements and f_surf in our model are mathematically analogous: the limiting current density of mass transport per catalyst area depends on the inverse of both the particle radius (r) and f_surf,^23,26,50 but neither of them affects the kinetic limitation.

As discussed in the previous section, the limiting current density of the V–T path depends on f_surf, thus the importance of the Heyrovsky step also depends on it. At sufficiently high values, the HOR is limited by the H₂ transport and the Heyrovsky step likely has little effect on the iη curve, whereas at small catalyst surface areas, its contribution could be significant. In the following, we analyze the effect of the Heyrovsky step in more detail using the dual pathway model that takes into account all three elementary steps.^1,15,28 The expressions of the kinetic model are given in the ESI of our earlier paper.¹ The concentrations, hydrogen coverage and overpotential of this reaction scheme were determined numerically as a function of i_el.

Fig. 7 shows the iη curves with low, but non-zero Heyrovsky exchange current density. The simulation parameters are the same as given in Fig. 6, except for r_H = 0.001 (r_H = ν_0,H/ν_0,V as defined by Wang et al.,¹⁵i.e. i_0,H = 0.1 mA cm⁻²). The aforementioned shoulder is clearly discernible in the iη curves, and occurs at current density i_el ≈ i_lim,HOR,VT (eqn (16b)) regardless of the value of f_surf (Fig. 7A). However, as Fig. 6B indicates, at extremely low catalyst surface ratios the shoulder may occur at such a low current density that noise in measured curves may limit their usefulness.⁴⁸ Up to the shoulder current density the V–T approximation is accurate, as the comparison of Fig. 6B and 7A shows, but it fails at current densities that exceed i_lim,HOR,VT.

Fig. 7 The HOR current densities normalized using (A) the total HOR limiting current density (eqn (17b)) and (B) the mass transport limitation simulation parameters are otherwise identical to those given in Fig. 6, except for r_H = 0.001.

At the lowest f_surf values (10⁻⁴ and 10⁻³) the V–T limiting current density (i_lim,HOR,VT) is achieved approximately at the same potential around 0.4 V vs. RHE. Because in the case of 10⁻⁴ the current density at this potential is approximately equal to the kinetic limiting current density of the V–T path (Fig. 3), the overpotential decreases with increasing f_surf. However, for all smaller f_surf values the overpotential is approximately constant, because they, too, are kinetics limited and the shoulder corresponds to the same current density per catalyst area. The constant shoulder overpotential might help to distinguish the kinetics limited measurements from those affected by the mass transport. The shoulder potential also naturally depends on the reaction kinetics and, for instance, in the results obtained by Chen and Kucernak⁴⁸ the shoulder occurred at about 0.1 V vs. RHE.

Fig. 8 shows how the catalyst surface fraction affects the dependence of the hydrogen coverage on the current density. The shoulder clearly appears here at i_el = i_lim,HOR,VT in all the cases, and as in the iη curves (Fig. 7), it is the smoothest at the lowest catalyst surface ratios, but turns sharply to a plateau at high f_surf values. The shape of the curve depends on the relative magnitude of the mass transport limit i_lim,H2 compared with i_lim,HOR,VT, similarly to the current density in Fig. 7. Increasing the catalyst surface fraction steepens the shoulder and improves the accuracy of eqn (16b) because a sufficiently high f_surf makes the shoulder current density mass transport limited, and prevents its further increase. In this case f_surf ≥ 0.1 are mass transport limited. Decreasing the Heyrovsky exchange current density would have a similar effect because the fraction of the total current density that corresponds to the V–T path would be increased, thereby improving the accuracy of the V–T approximation and eventually leading to results similar to those shown in Fig. 6.

Fig. 8 The effect of the catalyst surface fraction f_surf on the hydrogen coverage of the catalyst surface as a function of the normalized current density. The simulation parameters are the same as in the iη curves in Fig. 7. When f_surf ≤ 0.01, the mass transport limitation is higher than the total V–T limiting current density (the vertical line).

5. Conclusions

We showed the general applicability of 1D models for planar electrodes with a wide range of catalyst surface areas. Although our mass transport model does not describe very sparse catalyst arrays accurately, the inaccuracies are limited to the HOR near the limiting current density. At catalyst loadings that are practical for solar energy conversion, the 1D model is as accurate as the full 2D model. Phenomenologically, this means that the region of spherical mass transport near the particles can be neglected, which simplifies the mathematical analysis. In general, the operation of the planar nanoparticle array electrodes combines features from both the planar electrodes and the single-particle electrodes. We showed how these cases can be used in a complementary way to model and analyze the operation of the catalyst arrays by considering the current density per electrode or catalyst area.

Focusing on the special case of HER/HOR on Pt, we derived an analytical model for the electrode that also takes into account the surface hydrogen coverage of the catalyst. The usefulness of the hydrogen coverage as the starting point for the analysis of limiting current densities was also demonstrated. This method could be especially useful for reaction mechanisms that are limited by a combination of the reaction kinetics and mass transport similarly to HOR via the V–T path. Although this current density is not the limiting current density of the whole reaction when the Heyrovsky step also contributes to the current density, it can be observed as a shoulder or a plateau in the iη curve. By simplifying the expression of the current density of the Volmer step, we derived analytical expressions for the micropolarization and high overpotential ranges of the iη curve. Although neither of the simplified solutions describes the overpotential range from 50 mV to 100 mV well, both are accurate in the overpotential ranges, where their underlying assumptions about the overpotential are valid.

Because the 1D model is built based on analytical expressions, it is not only computationally light, but also readily applicable to the analysis of the experimental results. By showing the validity of the catalyst surface area-based approach for partially covered electrodes, we also showed that partially covered electrodes can be studied using 1D models for planar electrodes that couple mass transport and reaction kinetics together. In an example case, we demonstrated how the analysis is performed in practice by varying the exchange current density or the limiting current density of the mass transport, depending on the current density normalization used.

Acknowledgements

This work was funded by Nordic Energy Research within the project Nordic Initiative for Solar Fuel Development (NISFD).

References

1. E. Kemppainen, A. Bodin, B. Sebok, T. Pedersen, B. Seger, B. Mei, D. Bae, P. C. K. Vesborg, J. Halme, O. Hansen, P. D. Lund and I. Chorkendorff, Energy Environ. Sci., 2015, 8, 2991–2999.
2. E. Kemppainen, J. Halme, O. Hansen, B. Seger and P. D. Lund, Int. J. Hydrogen Energy, 2016 DOI:10.1016/j.ijhydene.2015.12.207.
3. M. Carmo, D. L. Fritz, J. Mergel and D. Stolten, Int. J. Hydrogen Energy, 2013, 38, 4901–4934.
4. P. C. K. Vesborg and T. F. Jaramillo, RSC Adv., 2012, 2, 7933–7947.
5. Fuel Cell Handbook, U.S. Department of Energy, Morgantown, West Virginia, 7th edn, 2004.
6. J. Lee, B. Jeong and J. D. Ocon, Curr. Appl. Phys., 2013, 13, 309–321.
7. J. R. Bolton, S. J. Strickler and J. S. Connolly, Nature, 1985, 316, 495–500.
8. S. Hu, C. Xiang, S. Haussener, A. D. Berger and N. S. Lewis, Energy Environ. Sci., 2013, 6, 2984–2993.
9. L. C. Seitz, Z. Chen, A. J. Forman, B. A. Pinaud, J. D. Benck and T. F. Jaramillo, ChemSusChem, 2014, 7, 1372–1385.
10. H. Döscher, J. Geisz, T. Deutsch and J. A. Turner, Energy Environ. Sci., 2014, 7, 2951–2956.
11. N. P. Dasgupta, C. Liu, S. Andrews, F. B. Prinz and P. Yang, J. Am. Chem. Soc., 2013, 135, 12932–12935.
12. Y. Chen, K. Sun, H. Audesirk, C. Xiang and N. Lewis, Energy Environ. Sci., 2015, 8, 1736–1747.
13. L. Ji, M. D. McDaniel, S. Wang, A. B. Posadas, X. Li, H. Huang, J. C. Lee, A. A. Demkov, A. J. Bard, J. G. Ekerdt and E. T. Yu, Nat. Nanotechnol., 2015, 10, 84–90.
14. H. A. Gasteiger, J. E. Panels and S. G. Yan, J. Power Sources, 2004, 127, 162–171.
15. J. X. Wang, T. E. Springer and R. R. Adzic, J. Electrochem. Soc., 2006, 153, A1732–A1740.
16. J. Durst, C. Simon and H. A. Gasteiger, J. Electrochem. Soc., 2015, 162, F190–F203.
17. M. Secanell, K. Karan, A. Suleman and N. Djilali, J. Electrochem. Soc., 2008, 155, B125–B133.
18. P. R. F. Barnes, A. Y. Anderson, J. R. Durrant and B. C. O’Regan, Phys. Chem. Chem. Phys., 2011, 13, 5798–5816.
19. W. Du, N. Xue, W. Shyy and J. R. R. A. Martins, J. Electrochem. Soc., 2014, 161, E3086–E3096.
20. T. J. Mills, F. Lin and S. W. Boettcher, Phys. Rev. Lett., 2014, 112, 148304.
21. K. Aoki, Electroanalysis, 1993, 5, 627–639.
22. H. Reller, E. Kirowa-Eisner and E. Gileadi, J. Electroanal. Chem. Interfacial Electrochem., 1982, 138, 65–77.
23. T. J. Davies and R. G. Compton, J. Electroanal. Chem., 2005, 585, 63–82.
24. M. A. Montero, M. R. Gennero de Chialvo and A. C. Chialvo, Electrochim. Acta, 2010, 56, 756–761.
25. A. Berger and J. Newman, J. Electrochem. Soc., 2014, 161, E3328–E3340.
26. A. J. Bard and L. R. Faulkner, Electrochemical Methods: Fundamentals and Applications, WILEY-VCH Verlag, 2001.
27. A. M. Chaparro, J. Electrochem. Soc., 2014, 161, E3078–E3085.
28. M. R. Gennero de Chialvo and A. C. Chialvo, Phys. Chem. Chem. Phys., 2004, 6, 4009–4017.
29. P. M. Quaino, M. R. Gennero de Chialvo and A. C. Chialvo, Phys. Chem. Chem. Phys., 2004, 6, 4450–4455.
30. P. M. Quaino, M. R. Gennero de Chialvo and A. C. Chialvo, Electrochim. Acta, 2007, 52, 7396–7403.
31. J. X. Wang, T. E. Springer, P. Liu, M. Shao and R. R. Adzic, J. Phys. Chem. C, 2007, 111, 12425–12433.
32. S. A. Vilekar, I. Fishtik and R. Datta, J. Electrochem. Soc., 2010, 157, B1040–B1050.
33. M. D. Arce, H. L. Bonazza and J. L. Fernández, Electrochim. Acta, 2013, 107, 248–260.
34. T. Shinagawa, A. T. Garcia-Esparza and K. Takanabe, Sci. Rep., 2015, 5, 13801.
35. M. R. Gennero de Chialvo and A. C. Chialvo, Electrochim. Acta, 1998, 44, 841–851.
36. K. Aoki, H. Toda, J. Yamamoto, J. Chen and T. Nishiumi, J. Electroanal. Chem., 2012, 668, 83–89.
37. P. J. Rheinländer, J. Herranz, J. Durst and H. A. Gasteiger, J. Electrochem. Soc., 2014, 161, F1448–F1457.
38. J. Durst, A. Siebel, C. Simon, F. Hasché, J. Herranz and H. A. Gasteiger, Energy Environ. Sci., 2014, 7, 2255–2260.
39. C. M. Zalitis, J. Sharman, E. Wright and A. R. Kucernak, Electrochim. Acta, 2015, 176, 763–776.
40. R. C. Rossi and N. S. Lewis, J. Phys. Chem. B, 2001, 105, 12303–12318.
41. C. Tsai, F. Abild-Pedersen and J. K. Nørskov, Nano Lett., 2014, 14, 1381–1387.
42. A. Bruix, J. A. Rodriguez, P. J. Ramírez, S. D. Senanayake, J. Evans, J. B. Park, D. Stacchiola, P. Liu, J. Hrbek and F. Illas, J. Am. Chem. Soc., 2012, 134, 8968–8974.
43. J. Zheng, Y. Yan and B. Xu, J. Electrochem. Soc., 2015, 162, F1470–F1481.
44. W. Sheng, H. A. Gasteiger and Y. Shao-Horn, J. Electrochem. Soc., 2010, 157, B1529–B1536.
45. CRC Handbook of Chemistry and Physics, ed. W. M. Haynes, CRC Press, 95th edn, 2014.
46. N. M. Markovic and P. N. J. Ross, Surf. Sci. Rep., 2002, 45, 117–229.
47. K. Elbert, J. Hu, Z. Ma, Y. Zhang, G. Chen, W. An, P. Liu, H. S. Isaacs, R. R. Adzic and J. X. Wang, ACS Catal., 2015, 5, 6764–6772.
48. S. Chen and A. R. Kucernak, J. Phys. Chem. B, 2004, 108, 13984–13994.
49. F. J. Vidal-Iglesias, J. Solla-Gullón, V. Montiel and A. Aldaz, Electrochem. Commun., 2012, 15, 42–45.
50. S. Chen and A. Kucernak, J. Phys. Chem. B, 2003, 107, 8392–8402.

---