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

Powered by TCPDF (www.tcpdf.org) This material is protected by copyright and other intellectual property rights, and duplication or sale of all or part of any of the repository collections is not permitted, except that material may be duplicated by you for your research use or educational purposes in electronic or print form. You must obtain permission for any other use. Electronic or print copies may not be offered, whether for sale or otherwise to anyone who is not an authorised user. Kemppainen, Erno; Halme, Janne; Lund, Peter


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,2The model system studied is a planar TiO 2 coated silicon photoelectrode (PE) partially covered with randomly distributed Pt nanoparticles (Fig. 1).4][5][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.
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.[13] For example, our previous work showed that a Pt loading as low as 100 ng cm À2 is enough to run the hydrogen evolution reaction (HER) in a PEC cell at 10 mA cm À2 current density without inflicting more than 50 mV overpotential. 1 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 À2 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 À2 . 3,5,14In this case, the catalyst surface area exceeds the projected surface area of the electrode by more than one order of magnitude. 15,168][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. 1 Our second paper extended the model by considering the transport of both gaseous and dissolved H 2 , as well as mass exchange between them in the electrolyte. 2In 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. 20his 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.

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 2 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 À2 . 2 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 driftdiffusion model, it can also be considered as a superposition of the overlapping, hemispherical flux profiles of the individual catalyst particles: [21][22][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. 23umerical 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,24Therefore, 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 This journal is © the Owner Societies 2016 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

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 ), 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 4 1, Fig. 2C), i.e. as the surface roughness factor. 15,25,26sing the ratio to describe surfaces that are only partially covered with the catalyst (f surf o 1, Fig. 2A) is apparently not as common, and we have found only a few examples of this. 24,27][30][31][32][33][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, 15 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.

Mass transport and surface concentrations
We assume linear mass transport to allow direct comparison with our earlier 2D simulations. 1 For the sake of simplicity and to allow analytical solutions to be reached, the effect of H 2 bubbles on the transport of the dissolved H 2 molecules is neglected here. 2n 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 (1 M) proton concentration and 1 bar H 2 pressure, i.e. 0.77 Â 10 À3 M, but it is possible that only a small fraction of H 2 is in a form that can participate in the reaction at the electrode surface. 36In our simulations, c b = c 0 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 mm diffusion layer i lim,H + = À35 940 mA cm À2 and i lim,H 2 = 15.4 mA cm À2 , see Table 1 for the simulation parameters).We follow the convention that positive currents correspond to HOR and negative currents to HER.
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, 27 we have neglected them here, as it seems that the HER/HOR on Pt could be described only by mass transport and reaction kinetics. 24

Reaction kinetics
The total hydrogen oxidation and evolution reaction proceeds through intermediate reactions called Volmer, Tafel and Heyrovsky steps H 2 $ 2H + + 2e À HER/HOR (3a) H 2 $ H ads + H + + e À Heyrovsky (3c) 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 39 as the dominant path.To allow direct comparison with our earlier 2D simulations 1 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 (y) as a function of the current density.Alternatively, we can also determine the current density as a function of y, as discussed in the next section.We begin with the reaction rates of the Tafel (n T ) and Volmer (n 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 b = 1/2 (and thus Z 0 /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 y.The exchange rates of the reactions are denoted by subscript 0 and the equilibrium potential of the reaction is V 0 (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 40 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,42qn ( 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 The exchange current density is i 0,V = FÁn 0,V , where F is the Faraday constant.When the current density and y 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 2 mass transport limitation, leading to a solution that depends only on the limiting current density of the H 2 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 2 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 2 transport. 16,43n 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 c b n c 0 n instead of c n c b n ), 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 2 concentrations, but this property and comparison to experimental results [37][38][39] are beyond the scope of our analysis.

Simulation parameters
Table 1 contains the simulation parameters introduced in the previous sections.We used the same values that we used previously 1 or the values derived from these simulations.The limiting current densities correspond to the solution of the 1D drift-diffusion problem over a 5 mm thick diffusion layer.In the case of protons, electroneutrality is enforced and perchlorate ions (ClO 4

À
) 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). 1 The radius of the simulation cell corresponded to a given catalyst mass loading (ng cm À2 , L cat ) The density of the catalyst material is r cat , which for Pt is r = 21 450 kg m À3 . 45Depending 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.

Current density and hydrogen coverage of the catalyst
As mentioned, eqn (5) includes three unknown variables, i el , y and Z 0 .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, y 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 = n 0,T /n 0,V (similarly to Wang et al. 15 ), the relative exchange rate of the Tafel step compared with the Volmer step.
Although one can consider r T 4 1 to correspond to the Volmer step as the RDS and r T o 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. 30ince 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, y 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 y 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 y A [0,1]) obtained using eqn (10) gives The numerator is the Tafel rate with the bulk H 2 concentration and the denominator corresponds to the effects of the H 2 transport.It is readily apparent that the proton concentration has no effect on the interdependence of the current density and y.
Using eqn (11) it can be shown that qi el /qy o 0 (when y 0 A (0,1) and y A (0,1)), so the extremes y = and y = 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 Although there are minor differences, this is equivalent to the solution of Wang and co-workers. 15,47The physically sensible solution corresponds to the negative root, because it yields y = y 0 , when i el = 0 (and c b H 2 ¼ c 0 H 2 ).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 c H 2 ¼ c 0 H 2 ), the highest current densities per catalyst area that the Tafel step allows are i lim;HER;kin;Pt ¼ À2r T i 0;V y 0 2 (13a) 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 surf i lim,HER,kin,Pt (14a) 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 4 1 and the term is even smaller than this fraction.
On the other hand, if r T o 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,H 2 -0, this expression becomes constant and depends only on the equilibrium coverage and the bulk H 2 concentration, which implies that if the bulk H 2 concentration is fixed by the experimental conditions, the H 2 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.

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 y, 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).
This journal is © the Owner Societies 2016 Because proton transport does not affect y 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), y -1 corresponds to i eli lim,HER,kin ) and the mass transport limit is not reached.In eqn (15) the kinetic limitation is neglected and therefore y -1 corresponds to i el -ÀN.
Because the mass transport losses are mainly dictated by the H 2 concentration (eqn ( 6) and ( 7)), the limiting current density of the H 2 transport is an important parameter, even when considering only the HER. 1,16,44However, 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 2 transport limitation.Moreover, because the H 2 transport and y 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 y -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,H 2 ) rewrite this as Setting y = 0 in eqn (12a), i.e.C = 0, yields the same result.(The other possibility that corresponds to C = 0, y = À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 2 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 2 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 26 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,49he 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 À4 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 E i lim,H 2 ), 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. 24 .

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 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, y is a function of the current density (eqn ( 12) and ( 15)) or current density is a function of y 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 y E y 0 and c H þ % c 0 H þ , 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 2 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.Z 0 /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 2 through the surface hydrogen coverage).

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 À2 (f surf E 5.3 Â 10 À6 ) to 1000 ng cm À2 (f surf E 0.53), except for the high HOR overpotentials and low catalyst 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 À4 ) to the mass transport control (blue, f surf = 10 À2 ), and the intermediate region where the current is limited both by the kinetics and mass transport (black, f surf = 10 À3 ).These three characteristic cases can be identified in the other figures as red, blue and black curves, respectively (Fig. 5-8).
This journal is © the Owner Societies 2016 loadings (1.0 ng cm À2 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 À2 (f surf E 5.3 Â 10 À4 ) 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.
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,26For a particle with a 5 nm diameter (D H 2 = 5.11 Â 10 À5 cm 2 s À1 and c H 2 = 0.77 mM) this limit is about 3.1 Â 10 4 mA cm À2 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 E 10 À3 (Fig. 3B).Considering that the HOR limiting current density of the 1D model becomes inaccurate between 1.0 ng cm À2 (f surf E 5.3 Â 10 À4 ) and 10 ng cm À2 (f surf E 5.3 Â 10 À3 , 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 24 ) 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.17)), with the full 1D solution of the current density-overpotential curve (iZ curve).In all cases, the micropolarization curve is accurate when |Z| o 50 mV, and the Tafel equations are accurate for |Z| 4 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 À4 -0.1.

HOR limiting current density
As eqn (7) illustrates, the H 2 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 2 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 o 10 À4 ) and high (f surf 4 10 À2 ) 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 E i lim,H 2 ) and significantly higher than the combined limit i lim,HOR,VT (Fig. 3).For example, in the case f surf = 10 À3 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 À3 (Fig. 6B), confirms that the HOR current density is limited by the combination of the reaction kinetics and mass transport.
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 iZ 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 iZ 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.

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,47nlike the Tafel rate, the Heyrovsky rate is directly affected by the applied potential; thus its rate can be increased with the overpotential. 15,28Therefore, 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 iZ curve. 15,28uch shoulders in the HOR current density curve have been observed experimentally at least by Chen and Kucernak 48 and Elbert et al.; 47 Wang et al. 15 interpreted the experimental results of Chen and Kucernak 48 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 2 transport and the Heyrovsky step likely has little effect on the iZ 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,28The expressions of the kinetic model are given in the ESI of our earlier paper. 1 The concentrations, hydrogen coverage and overpotential of this reaction scheme were determined numerically as a function of i el .Fig. 7 shows the iZ 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 = n 0,H /n 0,V as defined by Wang et al., 15 i.e. i 0,H = 0.1 mA cm À2 ).The aforementioned shoulder is clearly discernible in the iZ curves, and occurs at current density i el E 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. 48Up 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 .
At the lowest f surf values (10 À4 and 10 À3 ) 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 À4 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 48 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 iZ 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,H 2 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 Z 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.

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 iZ 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 iZ 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.
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 iZ curves in Fig. 7.When f surf r 0.01, the mass transport limitation is higher than the total V-T limiting current density (the vertical line).

Fig. 1
Fig. 1 Illustrations of real and simulated electrodes with 50 ng cm À2 Pt nanoparticle (d = 5 nm) loading having 2.7% surface area fraction (f surf E 0.027): (A) SEM image of Pt nanoparticles on TiO 2 .Reproduced from ref. 1 with permission from The Royal Society of Chemistry; (B) a scheme of the diffusion domain approach (r E 30 nm); (C) local current density distribution (mA cm À2 ) in the electrolyte for i el E À20 mA cm À2 .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.

Fig. 2
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.
This journal is © the Owner Societies 2016 Phys.Chem.Chem.Phys., 2016, 18, 13616--13628 | 13623 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

Fig. 5
Fig.5compares the simplified solutions for the micropolarization range and Tafel equation (eqn(17)), with the full 1D solution of the current density-overpotential curve (iZ curve).In all cases, the micropolarization curve is accurate when |Z| o 50 mV, and the Tafel equations are accurate for |Z| 4 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 À4 -0.1.

Fig. 5
Fig. 5 Comparison of the expression for the micropolarization range (m, 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 .

Table 1
The simulation parameters and their original source (if not based only on 1) Limiting current density of H 2 transport 15.4 mA cm À2 i lim,H + Limiting current density of H + transport À35 940 mA cm À2 c 0 in the calculations.Because no electrons cross the electrodeelectrolyte interface in the Tafel reaction, the Volmer rate multiplied by the Faraday constant and f surf yields the current density per electrode area. use This journal is © the Owner Societies 2016 Phys.Chem.Chem.Phys., 2016, 18, 13616--13628 | 13627