The Butler-Volmer Equation for Polymer Electrolyte Membrane Fuel Cell (PEMFC) Electrode Kinetics: A Critical Discussion

Most recent publications on PEMFC modelling describe the electrode kinetics of the catalytic reactions with a mechanistically simple equation usually called the “Butler-Volmer equation”. Here we argue that the Butler-Volmer equation is unsuited to PEMFC kinetics and offers no advantages over simpler empirical approaches. We also demonstrate that the most widely used formulation of this equation in the literature is inconsistent with the textbook equation. In order to critically assess the literature, we derive fundamental concepts underpinning the Butler-Volmer equation, without speciﬁc mechanistic assumptions. Through illustrative PEMFC simulations using various possible kinetic treatments, we demonstrate how models using an inappropriately formulated kinetic equation can give inaccurate predictions, or can contain superﬂuous detail. Since deviations in current distribution are often appreciableevenwhenthepolarizationcurvechangesnegligibly,experimentalvalidationofmodelsbycomparisontothepolarizationcurvealonecanbeinadequate.Forsimple,empiricalmodels,weadvocatelinearreversiblekineticsforthehydrogenoxidationreaction(HOR)andirreversiblekineticsfortheoxygenreductionreaction(ORR),avoidingsuperﬂuousparameterizationintheButler-Volmerequation.ReductionofempiricismrequiresmoremechanisticallydetailedmodelsthantheButler-Volmerequation.Weemphasizethatincorrectlyformulatedkineticequationsriskgeneratingincorrectdataormisleadingconclusions.©TheAuthor(s)2019.PublishedbyECS.ThisisanopenaccessarticledistributedunderthetermsoftheCreativeCommonsAttribution4.0License(CCBY,http://creativecommons.org/licenses/by/4.0/),whichpermitsunrestrictedreuseoftheworkinanymedium,providedtheoriginalworkisproperlycited.[DOI:10.1149/2.0361904jes]

The design and optimization of polymer electrolyte membrane fuel cells (PEMFCs) has historically been closely integrated with theory and simulation, in order to describe the many closely coupled physical phenomena occurring in different parts of the fuel cell device. 1 An especially intricate aspect of PEMFC modeling is the description of the catalyst layer, which traditionally consists of a composite of catalyst material (usually Pt supported on nanostructured carbon) and polymer electrolyte [2][3][4] and is the site of the hydrogen oxidation reaction (HOR) and oxygen reduction reaction (ORR) at the anode and cathode respectively.
During the last decade, significant developments have been made in the development of PEMFC physical models to account for mass and charge transport in systems with at least three phases (electrode, polymer electrolyte, and gas), as well as to describe non-isothermal and "two-phase" phenomena (incorporating liquid water). [5][6][7] Concurrently, several researchers have developed and advocated more detailed mechanistic models for the kinetics of the catalytic processes [8][9][10][11][12][13][14][15][16][17] as well as developing the innovative measurement techniques needed to interrogate electrode kinetics in a realistic context. [18][19][20] Nonetheless, as we shall demonstrate, the kinetic treatment in the majority of contemporary PEMFC models continues to use a mechanistically simple equation often referred to as the "Butler-Volmer equation".
In this paper, we provide a critical discussion of the widespread application of the latter class of 'simple' kinetic models in PEMFC simulation. Our key results are as follows: • We re-iterate the important but evidently often misunderstood point that the Butler-Volmer equation does not account, in general, for multistep phenomena in fuel cell catalytic reactions. Therefore, the Butler-Volmer equation is inherently an empirical treatment; although it has a basis in transition state theory for a one-electron process, 21,22 this rationalization does not apply in the fuel cell context. Previous authors have criticized the mechanistic basis for using the Butler-Volmer equation for fuel cell reactions, 23,24 and other authors acknowledge that its use constitutes a significant approximation. 25 Nonetheless, we will demonstrate through a review of contemporary literature that the limitations of the Butler-Volmer equation for multistep processes are still likely not universally understood.
• The "Butler-Volmer equation" in contemporary PEMFC literature often differs in its formulation from the textbook Butler-Volmer equation, in such a manner that the most commonly used variant of the equation fails to agree with the simple, empirical conception of total electrode current density being the sum of contributing oxidative and reductive current densities.
• Under common operating conditions, the electrode kinetics are insensitive to further addition of parameters in the Butler-Volmer equation, as compared to simpler linear kinetic relations, or the Tafel equation. Even when valid as an empirical treatment, the Butler-Volmer equation yields excessive parameterization. This leads to a valueless increase in the complexity of fuel cell models and hinders experimental validation.
• The term overpotential is widely misused with respect to its definition in classical textbooks and literature; we offer a simple, empirical derivation of the Butler-Volmer equation in order to express clearly the meaning of this quantity, and recommend alternative nomenclature for other related quantities often appearing in PEMFC kinetic equations.
• Polarization curve predictions are often 'blind' to incorrect formulation of electrode kinetics under typical PEMFC operating conditions, even when spatial current distribution still varies appreciably as a function of the kinetic model. We emphasize that fitting to an experimentally-derived polarization curve alone is often not sufficient to ensure the correctness of the choice of kinetic model.
• Conceptually doubtful kinetic models that appear suitable under conditions where the electrode environment is homogeneous can give erroneous predictions at higher current densities, or under extreme conditions where the overpotential varies in the plane of the electrode.
For these reasons, we argue that, in the PEMFC context, advantages due to the apparent generality of the Butler-Volmer equation are superficial -for circumstances where more detailed mechanistic models are deliberately eschewed, simpler approximations are preferred for both HOR and ORR.
We must stress that some of our determinations and recommendations have been made before, and we draw attention below to existing literature issuing similar warnings. It is likely that many PEMFC researchers are well aware of the above concerns and avoid the corresponding errors in their own work. In principle, much of the content of this work should be well known, or even obvious from standard physical chemistry principles; however, the systematic pattern of error and misunderstanding we have noted in contemporary publications suggests that prior remarks on electrode kinetics have not been heeded, and prompts us to consolidate these arguments in a single text focusing directly on PEMFC electrode kinetics.

What is the Butler-Volmer Equation?
By the 'textbook' Butler-Volmer equation -hereafter simply the "Butler-Volmer equation" -we refer to the common textbook form of an equation expressing the faradaic current density due to a redox process as the sum of independent contributions from oxidation and reduction. Both contributions include an exponential rate dependence upon electrode-electrolyte potential difference, whose symmetry between the oxidation and reduction directions is parameterized by empirical, unitless constants called transfer coefficients. Mathematically, by "Butler-Volmer equation" we mean: where E -electrode-electrolyte potential difference (V), named "electrode polarization" by Vetter 26 E eq -equilibrium potential, i.e. the value of E at which i = 0 f ≡ F/RT where F and R are respectively the Faraday and gas constants, and T is temperature (K) i -faradaic current density (A m −2 ) i 0 -exchange current density (A m −2 ) α,β -transfer coefficients (unitless) η -activation overpotential, or simply "overpotential" (V) We preface our future discussion by emphasizing that, in this expression, both exchange current density i 0 and equilibrium potential E eq depend on the local concentrations of reactants and products; this concentration dependence is implicit in Eq. 1, rather than being explicitly stated. We also emphasize the difference between the concentration-independent electrode polarization E and the concentration-dependent overpotential η, and their relation through Eq. 2.
From a historical perspective, it is curious to note that Erdey-Grúz and Volmer contributed to the development of this electrode kinetic equation through analysis of the hydrogen evolution reaction, 27 even though this reaction, like the HOR, is known not to always exhibit what today is called "Butler-Volmer kinetics". 24 Butler also worked on catalytic reactions of hydrogen, 28 though not initially. 29,30 In this work, we are not concerned with the historical development of the Butler-Volmer equation; rather, we reference its meaning as defined through standard textbooks of electrochemistry.
Suitable references giving a definition of "the Butler-Volmer equation" equivalent to Eq. Different choices of notation and mathematical expression between the above references might cause the mutual equivalence of their corresponding equations and their connection to Eq. 1 to be unclear. Their equivalence is discussed briefly below in Our Derivation of the Butler-Volmer Equation and Comparison with Bernardi-Verbrugge Formulation. A full discussion of the interrelation between different forms of the Butler-Volmer equation and the advantages and disadvantages for selection among them is an intended topic for our future work, but has relevance also outside the PEMFC field and exceeds the scope of this article.

The Need for an Accurate Kinetic Equation
To express the local performance of the catalyst layer (spatially resolved, but averaged over length scales larger than microstructural particle sizes) in a continuum mathematical model of a fuel cell, the minimal requirement for the electrode kinetic equation is that it should yield the volumetric current density due to the faradaic reaction in the catalyst layer (A m −3 , per volume of catalyst layer) as a function of the electrode-electrolyte potential difference and the concentration(s) of the reactant(s), and, for a reversible process, also the concentration(s) of the product(s). In a non-isothermal model, temperature variation must also be considered.
Crucially, the choice of kinetic equation will not only impact the kinetically-controlled regime of PEMFC operation that occurs at low activation overpotential where the current density is small. Because the kinetic equation links the current density to the reactant concentration, the kinetic equation is also needed to predict the contribution of any transport-controlled limitation. This can be explained as follows: if the kinetic equation lacked any concentration dependence, no amount of depletion of reactant would alter the predicted current density. Rather, in a high current density limit, a concentration-independent kinetic equation would predict physically absurd results where the current density will be non-zero even in the presence of zero or 'negative' concentrations of reactant.
Hence, any error in the concentration dependence of the chosen kinetic equation may yield incorrect predictions at higher current densities, even if the predictions are accurate at low current densities where the reactant concentration is approximately uniform. The importance of the kinetic equation for the accuracy of predictions in the transport-controlled regime was emphasized by Moore et al. 31 In general, spatial non-uniformity of current distribution (and hence electrode utilization) in a fuel cell implies a non-uniform overpotential in the plane of the electrode. Simulations that yield data on spatial non-uniformity depend on the kinetic equation to accurately predict the overpotential; if the electrode kinetic equation is incorrectly formulated, it is likely that predictions concerning spatial inhomogeneities will also be unreliable.

Standard Approaches in PEMFC Electrode Kinetics: An Overview of Historical and Contemporary Literature
In order to make an informed judgement on trends in the description of electrode kinetics in the literature on modeling for PEMFCs, we conducted a comprehensive literature review, covering 276 PEMFC modeling papers from 1989 to 2018. Within these, we also focus on the subset of 34 papers published from 2017 to 2018, in order to compare contemporary practice to overall historical context.
The review considers only spatially resolved models in which a variation of reactant concentration is allowed in the PEMFC: either with spatial distribution, or in the simplest case through a predicted transport loss between inlets and the site of reaction. We have chosen to limit the scope of our review to fuel cell studies considering conventional catalysts and materials: Pt combined with a carbon-based conductive additive and Nafion ionomer, at temperatures typically close to T = 80°C. Consequently we have ignored papers relating to hightemperature PEMFCs and novel membrane or catalyst materials. We have also focused principally on research papers predicting polarization curves and/or losses in the catalyst layer, or otherwise modeling the performance of the fuel cell as a device, rather than simply constructing a microscopic model of one feature. In general, we have included only papers citing a full equation set, and have not included those that directly duplicate a previously reported model for a new geometry, as in a number of flow channel configuration studies. We have also have excluded those papers whose mathematical description is insufficiently complete to classify according to the kinetic relations given below, except where missing information can be inferred safely by comparison to related cited or citing articles.
Given the extent of the review, it is very likely that we have missed some publications, and that there exist some occasional errors both in our transcription and interpretation as well as in our rigour in applying the above limitations; we believe that any such errors do not affect our overall judgement on trends. Comprehensive details of the review,  The review identified three common classes of electrode kinetic model, defined as follows: 1. Springer model -a Tafel law for the ORR (Eq. 3, below), combined with either no model for the HOR (cathode-only simulation) or an assumption that the HOR is infinitely fast and so the Nernst equation (η = 0) can be used. This approach is so-named since it is common among models drawing primary influence from the foundational 1991 work by Springer et al. 32 2. Wang-Um model -a Tafel law for the ORR, combined with a linear equation in electrode polarization for the current density due to the HOR (Eq. 4, below). This approach has become standard among models influenced by the 2000 work of Um et al., 33 and appears frequently in other subsequent works from the group of C.Y. Wang. 34 3. Bernardi-Verbrugge formulation -an adapted Butler-Volmer equation (Eq. 5, below) at one or both electrodes; the nonequivalence of this equation with the Butler-Volmer equation will be discussed further in the following section. This particular and widely repeated equation form originates in the major 1991 and 1992 works of Bernardi and Verbrugge, hence our name. 35,36 We avoid the term "model" in this case since, although the equation differs from the standard Butler-Volmer equation, we believe that the original authors did not intend to construct a distinct physical model. a [5] Quantities denoted "ref" are constants. In the Bernardi-Verbrugge formulation Eq. 5, m denotes either anode or cathode, with c m = c H 2 at the anode and c m = c O 2 at the cathode.
Miscellaneous models not fitting the above classification include 'true' Butler-Volmer models, models invoking detailed (multistep) chemical mechanisms, and a few rare cases such as the use of the irreversible Tafel law to describe the HOR.
Among the above models, the Springer and Wang-Um models are both simple descriptions, involving clear approximations of slow, irreversible kinetics for the ORR, and fast, reversible kinetics for the HOR. Accordingly, they also require a minimal number of empirical parameters to express the kinetic dependences upon polarization and reactant concentrations. In contrast, the Bernardi-Verbrugge equation, as well as more mechanistically detailed models, require more extensive parameterization. Figure 1 illustrates the historical distribution of kinetic models among the classes defined above. The Springer and Wang-Um models combine to give 42% of modeling works; the Bernardi-Verbrugge formulation also accounts for 42%. Models using a 'true' Butler-Volmer equation or a detailed multistep mechanism are scarce in general; apart from works considering contaminants, in the whole history of PEMFC modeling we identified only twelve full cell PEMFC simulation papers including a multistep kinetic model at one or both electrodes; 14  work from the groups of M. Secanell and A.Z. Weber accounts for more than half of these. Figure 2 focuses specifically on publications dated from 2017 and 2018. In this sample of the contemporary literature, the Bernardi-Verbrugge formulation is used in the majority (65%), while simpler modeling approaches have become less frequent. One possible cause for this development is that the Bernardi-Verbrugge formulation has been implemented as a built-in equation in the increasingly popular commercial simulation software ANSYS Fluent, 48 49 and also implies that a principal inspiration is the work by Mazumder et al., 51 which uses the Bernardi-Verbrugge formulation.
In the literature, the Bernardi-Verbrugge formulation Eq. 5 is often called the "Butler-Volmer equation" (in dozens of papers, but see for  52 ); however, it is not equivalent to the Butler-Volmer equation Eq. 1. This is because Eq. 5 contains only one concentration-dependent term, in the common multiple preceding the two-term bracket. This confusing feature means that the concentration dependence applies equally to both forward and reverse (oxidative and reductive) processes. To our knowledge, this is not compatible with any known mechanistic conception of either the HOR or ORR, and was highlighted for criticism in the insightful but apparently largely overlooked 2009 work by Andreaus et al., who wrote: "Unfortunately, in the literature on fuel cell modeling, concentration dependencies are often written… as if they would apply equally to forward and back reactions." 3 Also, in the limit of i ref tending to infinity, the Bernardi-Verbrugge formulation does not yield the Nernst equation relating equilibrium potential to reactant and product concentrations, as occurs in the corresponding limit of the Butler-Volmer equation and should be expected for a kinetic model that is thermodynamically consistent with an infinitely fast, reversible process.
We emphasize that in raising this issue with the Bernardi-Verbrugge formulation we do not wish to single out these two authors for criticism; the original studies 35,36 marked a major development in considering fuel cell losses arising from a wide range of physical phenomena. Our concerns apply to the later, uncritical re-use of these equations, in a range of regimes.
The problems in the Bernardi-Verbrugge formulation can be clarified by direct and more detailed comparison with the Butler-Volmer equation. To facilitate this, we present here -as briefly as possible -an empirically motivated derivation of the Butler-Volmer equation based on a minimal set of assumptions about the electrode kinetic phenomena.

Our Derivation of the Butler-Volmer Equation and Comparison with Bernardi-Verbrugge Formulation
To derive the Butler-Volmer equation, we make the following assumptions: 1. The total current density due to a reversible electrochemical reaction can be expressed as the sum of an oxidative and reductive contribution, which can be specified independently. 2. Both the oxidative and reductive contributions obey a Tafel equation, such that their rates are directly proportional to an exponential function of the electrode polarization; that is, for constant reactant and product concentrations maintained by infinitely fast mass transport, they exhibit some empirical Tafel slope with respect to polarization of the electrode. Empirical Tafel slopes can be expressed as a transfer coefficient, according to IUPAC guidance. 53 The Tafel slopes are not assumed to be correlated between the two reaction directions. 3. Each of the oxidative and reductive contributions is proportional to its own dimensionless multiplying factor g j , which expresses the concentration dependence of the reaction. The quantity g j is defined to equal 1 under some set of reference concentrations, and tends to zero as the reactant concentrations go toward zero. 4. At zero net current density, a condition of dynamic equilibrium between anodic and cathodic processes exists, such that the polarization under these conditions (that is, the equilibrium potential) obeys the thermodynamic Nernst equation.
Note that we do not make any assumptions regarding whether the process is a multistep process or whether it is an elementary step with a single electron transfer. Although the exponential dependence of current density upon polarization in the Tafel equation is often justified on the basis of transition state theory in the case of an elementary step, 21,22 neither the HOR nor the ORR is an elementary step. Hence for the multistep processes of the HOR and ORR, the weakest assumption is Assumption 2 (a rate law given empirically by the Tafel equation), since it implies that an electron transfer step is always rate-limiting; in a multistep process combining electron transfer with other kinetically limited chemical processes such as adsorption, there is no reason to expect this to necessarily be the case. We do not wish to argue that this assumption is true for PEMFC electrode kinetics; rather, we need to consider the consequences when it holds in order to assess the validity of the Butler-Volmer equation and its related Bernardi-Verbrugge formulation. Unlike the derivations of the one-electron Butler-Volmer equations in standard textbooks, 21,22 we do not invoke any a priori assumptions about the electrochemical mechanism -our derivation is empirical.
Assumptions 1-4 can be expressed mathematically for a reversible electrochemical process m as follows: All quantities are physical constants of the reaction, except for polarization E and the concentration dependences We denote this equilibrium potential as E eq,ref,m , so that, since all g j,m = 1: The quantity i 0,ref,m defined by Eq. 9 is the exchange current density under the reference conditions, where by "exchange current density" we mean the implied absolute contribution from oxidative or reductive current density under equilibrium conditions. By substitution of Eq. 9 into Eq. 7 and Eq. 8 and then combining using Eq. 6: [10] For brevity of notation, we introduce a concentration-independent polarization with respect to reference conditions as: Eq. 12 has the recognisable form of the Butler- At present we have introduced no correlation between the Tafel slopes or the empirical concentration dependences for the two reaction directions. Such a correlation is required, however, when considering deviation from the reference conditions. Under any conditions distinct from the reference conditions -that is, for any defined concentrations c i -there exists an equilibrium condition where E = E eq,m and i m = 0. So, E eq is a general function of c i , as known from the Nernst equation. From Eq. 10, we write for i m = 0: Rearranging: g ox,m g red,m exp (β m + α m ) f E eq,m − E eq,ref,m = 1 [14] or The thermodynamic Nernst equation requires, however, that: ln g eq,m [16] where g eq,m is a certain function of the normalized concentrations, depending in the ideal case only on the stoichiometry of the reversible reaction when written as a one-electron transfer, but in general expressing activity effects also. We define g eq,m to equal 1 under the reference conditions. Subject to our assumptions, both Eq. 15 and Eq. 16 must hold simultaneously -hence, there exists a certain constraint in order for the initially assumed Tafel laws to be thermodynamically compatible with a reversible process: g ox,m g red,m = g eq,m (βm+αm ) [17] That is, if the Tafel slopes are known at constant concentration, the microscopic reversibility of any reaction described by a Butler-Volmer equation requires a certain relationship between the empirical concentration dependences of the oxidative and reductive currents: the apparent reaction orders cannot be specified independently while retaining thermodynamic consistency. It is useful to recognize that no such constraint applies to a reaction considered irreversible, since the thermodynamic equilibrium of such a reaction is unattainable. An extremely common special case of Eq. 17 is that of an elementary step where a single molecule A undergoes first-order reversible one-electron reduction to a single molecule B. If activity effects can be ignored, then in this case g eq = g ox / g red , and Eq. 17 reduces to β m = (1 -α m ); this substitution is already included in the forms of Eq. 1 defined within the referenced standard textbooks. 21,22,26 Depending on whether the oxidative or reductive concentration dependence is specified a priori, two possible forms of the Butler-Volmer equation arise from combination of Eq. 17 with Eq. 12: Comparing Eq. 18 or Eq. 19 with the Bernardi-Verbrugge formulation Eq. 5, it is clear that the latter treatment of PEMFC electrode kinetics implies that g eq,m = 1, within the interpretation of total current density being the sum of oxidative and reductive contributions. This is only true if there is negligible deviation from reference conditions; however, it is very rare for this approximation to be made explicit. Insofar as the deviation from reference conditions is assumed negligible, any expression of concentration dependence in the Bernardi-Verbrugge formulation is also incomplete -as noted in The Need for an Accurate Kinetic Equation section above, this would make the description unsafe in general under transport-controlled conditions. A similar observation was made by Andreaus et al. in 2009: "…the original Butler-Volmer equation … should be used only when electrode potential and all concentrations are uniform. Such conditions are barely encountered in fuel cell electrodes." 3 The most general definition of the Butler-Volmer equation arises from further consideration of the equilibrium potential as expressed in the Nernst equation. Using Eq. 16, the electrode polarization can be written as follows: For any concentrations c i , we can write the overpotential -that is, the polarization with respect to the equilibrium potential at which i m = 0 -as: [21] which therefore gives, from Eq. 20: ln g eq,m [22] Using Eq. 22, we can expand the exponential terms of the Butler-Volmer equation as: +αm [24] to give, from Eq. 18: The concentration dependence is now shared by both terms in the sum and can be factorized out, to give a general Butler-Volmer equation as: where the general exchange current density is: Eq. 26 is the same equation as the one defined by us as the Butler-Volmer equation at Eq. 1 and also derived through a variety of means in standard textbooks. 21,22,26 Here, we have shown it is a necessary consequence of the simple, empirical assumptions at the beginning of this section. Crucially, the Butler-Volmer equation, Eq. 26, contains two concentration-dependent terms: the exchange current density and the overpotential. Likewise, the Butler-Volmer equation in the "electroanalytical" form Eq. 12 contains two concentration-dependent terms: g ox,m and g red,m . The appearance of two independent concentration dependences reflects the two independent contributions to the current density: the forward and reverse reactions. The Bernardi-Verbrugge formulation Eq. 5, which has just one concentration dependence, cannot reproduce the same behavior.
The simplest interpretation of the Bernardi-Verbrugge formulation is that it is an approximation to the Butler-Volmer equation in which η ref,m has been used where η m is required: that is, the concentrationindependent polarization has been used in place of the concentrationdependent overpotential. This is equivalent to neglecting deviations in the equilibrium potential due to the variation of local concentrations from their reference values -the approximation implied by the formulation is erroneous where such variations become appreciable.

Why Do "Incorrect" Kinetic Equations Give Valid Predictions?
The most commonly used kinetic description in contemporary PEMFC models is the Bernardi-Verbrugge formulation. We have shown that it is erroneous to call this formulation "the Butler-Volmer equation", because it does not equate to the textbook Butler-Volmer equation. We have further shown that the latter is a requirement of a simple and rational empirical statement considering the concentration dependence and thermodynamic consistency of the electrode kinetics; therefore, the Bernardi-Verbrugge formulation is incompatible with such statements. However, the Bernardi-Verbrugge formulation has, in diverse publications, yielded results that successfully predict the performance of a fuel cell as measured from the polarization curve, at least within certain regimes of operation (see, for example, Refs. 36,[54][55][56][57]. In this sense the Bernardi-Verbrugge formulation certainly seems to be valid for PEMFC electrode kinetics.
These observations can be rationalized by considering the typical kinetic performance of the two catalytic processes. Under typical fuel cell operating conditions, the total activation loss is dominated by the cathode. That is, relatively speaking, the anode kinetics are very fast, while the cathode kinetics are very slow.
For a very slow process (such as the ORR) under conditions where appreciable current density is drawn, it is suitable to consider the forward reaction as irreversible, so that the oxidative contribution to the reversible kinetic equation for the ORR is approximately zero and the activation overpotential is large. In these conditions, the Bernardi-Verbrugge formulation gives: [28] while the Butler-Volmer equation gives a Tafel law (directly from Eq. 12): In consideration of which terms are concentration-dependent and which terms are constant, these equations are identical. The Bernardi-Verbrugge equation therefore describes the cathode kinetics accurately because the part of the equation which disagrees with the Butler-Volmer equation is approximately zero under the most common conditions of operation. Of course, in this case, the generality of including the backward reaction is superfluous, and applying the irreversible Tafel equation (Eq. 29) directly, as in the Springer and Wang-Um models, is simpler and equally descriptive. For compatibility between Eqs. 28 and 29 it is also essential that the pre-exponential reaction order γ O 2 matches the empirical concentration dependence of g red,c and not the concentration dependence in the general expression for the Butler-Volmer exchange current density given by Eq. 27; in fact, this condition was not heeded in the original Bernardi-Verbrugge work, and the resulting 'hyper-correction' to an already correct Tafel law may account for the use of reaction orders close to 0.75 in some models. 58 Therefore, it is incorrect to speak of the pre-exponential term in the Bernardi-Verbrugge formulation as the "exchange current density". A clear summary of this concern was given in the remark by Djilali et al. in 2008 that, "The form of the Butler-Volmer equation used in current PEMFC models corresponds to forward reaction dominated processes and needs to be revisited in the case of backward dominated reactions." 59 In contrast to the Tafel law limit, for a very fast process (such as the HOR), we can consider the limit of high i ref,a . Then, the exponential terms in the Bernardi-Verbrugge formulation can be linearized to: [30] which is the linear form used in the Wang-Um model, while the Butler-Volmer equation gives through a similar linearization at high i 0,a : These equations are not exactly equivalent, due to the presence of the g eq,a terms in the latter. Under the approximation of negligible deviation from reference conditions, however, we can set g eq,a ≈ 1, and we have from both Wang-Um and Bernardi-Verbrugge: [32] while from Butler-Volmer: [33] So, the two descriptions become equivalent in their prediction of polarization, but only in the case where the anode H 2 concentration is approximately uniform. Where this requirement does not hold and the H 2 concentration becomes appreciably non-uniform for any reason, neither the Wang-Um model nor the Bernardi-Verbrugge formulation could be expected to give an accurate prediction of the spatial profile of either H 2 concentration or anodic current density in the anode catalyst layer. In both cases, the problem is associated with the use of polarization where overpotential is required.
We note that in the Wang-Um model's definition, the use of polarization in place of overpotential is not invoked in the summary given in the review paper by Wang et al. 34 but is invoked in the actual model defined by Um et al., 33 as well as all subsequent models we consulted. Um et al. is also cited by Li et al. for their concentration-independent definition of overpotential in the Bernardi-Verbrugge formulation used within the original Fluent model, 49 although Um et al. did not advocate this substitution for a nonlinear, reversible equation. We will use "Wang-Um model" to refer to anode kinetics linear in polarization (defined by Eq. 4) as commonly used by models in the literature, and use the terminology "corrected Wang-Um model" to refer to models with a linearized Butler-Volmer equation as given by Eq. 31.
Because the Wang-Um model and Bernardi-Verbrugge formulation equate to the Butler-Volmer equation for anode kinetics only when the H 2 concentration is near uniform, the concentration dependence included in either of the former descriptions lacks any physical meaning. Inconsistencies are expected whenever spatial variation of H 2 concentration becomes appreciable. Inconsistency might also be expected in the prediction of anode polarization as a function of position, but the usual high value of i ref,a would make these discrepancies small on an absolute scale.

Probing Possible Failure Regimes of Common Kinetic Treatments
From the analysis above, the Bernardi-Verbrugge formulation of cathode electrode kinetics is expected to give similar results to the true Butler-Volmer equation only if the overpotential is large, and the pre-exponential concentration dependence is identical for both formulations. In these standard operating conditions, we expect that the corrected Wang-Um model will be equally precise as the true Butler-Volmer equation, while requiring fewer parameters.
Also, significant variation in H 2 concentration at the anode reaction sites -as might occur in anode starvation conditions, for example -will cause the Bernardi-Verbrugge and Wang-Um descriptions to deviate from the predictions of the Butler-Volmer equation.
To investigate these cases, we performed two suites of simulations: a) a 1D model to confirm our assessment of the likely equivalence between the Bernardi-Verbrugge and Butler-Volmer equations when the cathode overpotential is large, and their mutual equivalence to the corrected Wang-Um model under standard operating conditions; b) a 2D model to scrutinize discrepancies between Wang-Um, Bernardi-Verbrugge and Butler-Volmer equations arising under anode starvation conditions. A relatively simple fuel cell model was employed, in order to emphasize the discrepancies due to the choice of kinetic model.
Full details of the model are given in Supplementary Material 3; the key assumptions are: Model (a): • 1D geometry parallel to the direction of the proton current through the membrane, including only the catalyst layers and membrane.
• Mass transfer of gases according to Fick's law.

Model (b):
• 2D geometry in a plane normal to the bulk gas flow, including gas diffusion layers, catalyst layers, and membrane.
• Momentum transfer for gases in porous regions according to the Darcy-Brinkman equations.
• Mass transfer of gases according to the Maxwell-Stefan equations, with a Bruggeman correction for tortuosity applying to all Maxwell-Stefan diffusivities.
Common to both models: • Charge transfer by electron conduction in the electrode material according to Ohm's law. • Charge transfer by proton conduction in the polymer electrolyte membrane according to Ohm's law.
• Kinetic equations applied directly with respect to gas concentrations in the catalyst layer, assuming thermodynamically controlled dissolution of the gases into the polymer electrolyte, according to Henry's law. No agglomerate model included.
• System assumed to be isothermal.
• Two-phase phenomena ignored, including electroosmotic drag through the membrane.
We acknowledge that the use of Ohm's law in the membrane, and the lack of consideration of any electroosmotic drag or two-phase phenomena, are incompatible with a quantitative prediction of fuel cell performance -as noted above, our goal is not a precise model but rather to emphasize the impact of the kinetic model and its interaction with conditions of non-uniform concentration of reactant in the catalyst layers. Discrepancies in spatial distribution of current density and reactant concentration due to the interaction of the kinetic equation with mass transport apply a fortiori to two-phase, non-isothermal models with more detailed and more realistic descriptions of the polymer electrolyte membrane and catalyst layers. Insofar as this work seeks to challenge inadequately appraised or superfluously complex equations in PEMFC models, we do not wish to cloud the description of kinetic phenomena by including unnecessary detail elsewhere in the mathematical description.
At the anode in the 1D model, linearized Butler-Volmer kinetics are used; in the 2D model, a range of electrode kinetic descriptions are compared. All models are based on transfer coefficients α a = 1 2 , β a = 1 2 and concentration dependences g ox,a = (c H 2 /c H 2 ,ref ) 1 2 , g red,a = 1. These settings retain the apparent concentration dependence of (c H 2 /c H 2 ,ref ) 1 2 at higher overpotential, as in many implementations of the Bernardi-Verbrugge formulation as well as the Wang-Um model (see, for example, 60,61 ); however, unlike the value of α a + β a = 2 in most implementations of the Wang-Um model, 34,60 the transfer coefficients used here are consistent with the constraint Eq. 17 for a thermodynamically compatible Butler-Volmer equation, so long as activity effects are ignored in the Nernst equation and so g eq = (c H 2 /c H 2 ,ref ) 1

.
Since the Wang-Um model gives a linear electrode kinetic equation, the change in transfer coefficient amounts simply to a rescaling of the magnitude of the reference current density at the anode. In the comparative implementations of the Wang-Um model and Bernardi-Verbrugge formulation, the empirical reaction order γ H2 is set to 1 2 . For the true Butler-Volmer equation, we set g red,c = (c O 2 /c O 2 ,ref ) (first-order kinetics in oxygen concentration), g ox,c = 1 (assumed uniform availability of protons and water), α c = 1, β c = 3. The latter transfer coefficient is determined by the constraint Eq. 17. The reference current densities are given a typical large value for the HOR (0.01 A cm −2 ) and a small value for the ORR (2.5 × 10 −6 A cm −2 ).
For the Butler-Volmer equation, we implemented the current density condition both in the overpotential form Eq. 26 and in the electroanalytical form Eq. 12. For the specific cathodic reaction behavior in question these give, respectively: [35] and so [37] Eqs. 34 and 35 are mathematically equivalent, so the results should be identical, although one form or the other might have better computational performance and reliability. In the 2D model where water vapor partial pressure is resolved locally in the catalyst layer, Eq. 35 is adapted to include a concentration dependence upon water for the oxidative contribution according to the constraint Eq. 17. However, since this term is always near zero, it makes a negligible change to the results. For the Bernardi-Verbrugge formulation, we considered the case Eq. 38 where the pre-exponential concentration dependence follows the Tafel law reaction order (first-order). At high activation overpotential, this treatment is only in error in the negligible oxidative contribution.
We also considered the same Bernardi-Verbrugge formulation when the O 2 reaction order is erroneously 'corrected' from 1 to 0.75 by analogy to Eq. 34, even though the Bernardi-Verbrugge formulation uses polarization rather than overpotential and so the resulting equation now differs mathematically from the others. Figure 3 shows the predicted polarization curves from the 1D model (a), comparing the Butler-Volmer equation to the Bernardi-Verbrugge formulation. The two mathematically equivalent Butler-Volmer equations and the inequivalent Bernardi-Verbrugge formulation (γ O 2 = 1) all give current density values that are indistinguishable within the numerical tolerance, which corresponds to equality to within five or more significant figures for the current densities achieved at E cell < 0.95 V. For the erroneous γ O2 = 0.75, a small discrepancy can be observed in the polarization curve at higher current densities, where the O 2 concentration varies to a greater extent. This deviation is well within the range we judge as typically acceptable for an experimentally valid fit, but might be corrected for fitting purposes by modification of one or more other physical parameters of the fuel cell -of course, such a 'correction' to the fit is needed only due to the incorrect kinetic equation, which therefore yields an incorrect value for one or more physical quantities in the model, whenever the latter are fitted to the polarization curve.  Figure 4 compares the polarization curves between the cases of a Butler-Volmer equation and a Tafel equation at the cathode -the latter model is equivalent to the corrected Wang-Um model according to our classification above. In the kinetically-controlled regime, the current densities are higher for the Tafel equation due to the absence of a reverse current contribution, as shown in the inset to Figure 4; again, however, the two equations agree to five significant figures at E cell < 0.95 V, because the cathodic contribution dominates completely. This makes clear that the equivalence of the Bernardi-Verbrugge and Butler-Volmer descriptions in Figure 3 arises because both are indistinguishable from the simple Tafel law, as predicted in the analysis in the section Why Do "Incorrect" Kinetic Equations Give Valid Predictions? above. Figure 4 confirms that the full Butler-Volmer equation, even if implemented correctly according to its textbook definition, introduces one superfluous parameter (an additional transfer coefficient) without adding any accuracy to the predictions of the model under normal conditions of operation. Although the results differ at very low current density, the mechanistic complexity of the real ORR compared to the simple assumptions implied by the Butler-Volmer equation means that the difference in prediction from the more complex Butler-Volmer equation does not imply any improvement in accuracy over the Tafel law. Figure 5 compares the predicted concentration profiles of O 2 for the 1D model (a) according to the different kinetic models. The concentration profile is plotted along a line through the thickness of the catalyst layer, in the 10 μm thickness closest to the boundary with the gas diffusion layer (GDL) where there exists a depletion zone for O 2 in the transport-controlled regime (data at E cell = 0.5 V). These data again demonstrate that while the Bernardi-Verbrugge and Butler-Volmer models are indistinguishable when they both reduce to the same Tafel equation, there is a noticeable difference in the spatial predictions of the Bernardi-Verbrugge formulations in the case of an erroneously 'corrected' pre-exponential reaction order. The difference in the concentration profile is much more marked than the difference in the polarization curve, underlining that the polarization curve is often not a reliable metric for validation of a theoretical model, when considered in isolation.
In the 2D model (b) (see Supplementary Material 3 for a diagram of the geometry), we compared the Wang-Um, Bernardi-Verbrugge, and Butler-Volmer descriptions for a model in which the anode is fed with fully humidified 20:80 H 2 :N 2 , rather than fully humidified pure H 2 as in the 1D model. We expect that these non-standard substoichiometric conditions will provoke transport-limiting behavior at the anode, rather than at the cathode as in typical PEMFC operation.  Figure 6 shows the polarization curves from the three cases. All are very similar, due to the small activation overpotential for the HOR; overall, the Butler-Volmer equation tends to differ from the Wang-Um and Bernardi-Verbrugge equations for intermediate current densities, since both of the latter use polarization in place of overpotential. Small differences arise between all three cases for E cell < 0.75 V (inset to Figure 6) where the H 2 partial pressure falls close to zero at the outlet end of the catalyst layer, creating a progressive redistribution of the anodic reaction and initiating more strongly transport-controlled behavior.
While the polarization curves do not show a marked difference between the kinetic models, the spatial distributions of reactant concentrations and anodic current density vary much more significantly. Figure 7 plots the concentration profile of dissolved H 2 in the polymer electrolyte phase through the thickness of the anode CL, in a region where H 2 depletion is expected, 0.2 mm upstream from the outflow.   Figures 8 and 9 plot the local volumetric current density within the anode CL: Figure 8 considers the variation through the CL thickness, as for the concentration profile; Figure 9 considers the variation in the cross-channel plane along the GDL-CL boundary from the anode inlet to outlet face (that is, along the line x = w GDL as defined in Supplementary Material 3). For both the concentration and local reaction rate data, it is apparent that the three different kinetic models applied to the anode give substantially different predictions -the distributions are, in some cases, qualitatively altered.
These results illustrate that an inappropriate kinetic description is likely to introduce unreliable spatial predictions whenever spatial variation of concentrations becomes important, even when the polarization curve is near-identical between different kinetic models. Spatially varying conditions are especially relevant in non-standard circumstances of particular concern for PEMFC performance and safety, such as high current densities, anode starvation, and transient response including cold start. Since spatial distribution of current density (and  hence catalyst utilization) controls spatial distribution of heat source and water source/sink, such considerations are perhaps even more significant in more physically complex models such as the two-phase, non-isothermal PEMFC models common in current practice.

Recommendations for Future Modeling Work
It is doubtful that the 'true' Butler-Volmer equation is mechanistically precise for PEMFC catalytic processes, given the general incompatibility of the multistep nature of the HOR and ORR with the approximations used in its derivation above. We also note that the apparently general Bernardi-Verbrugge formulation in fact equates to the empirical Butler-Volmer equation derived above in Our Derivation of the Butler-Volmer Equation and Comparison with Bernardi-Verbrugge Formulation only under the conditions of: a) negligible overpotential at the anode; b) high overpotential at the cathode. In these cases, the presence of independent forward and backward reaction terms in the electrode kinetic equation is unnecessary, and a simpler empirical approach can be adopted with a reduced requirement for parameterization.
It is also possible that the unphysical predictions of the Bernardi-Verbrugge formulation may lead to greater numerical instability in circumstances where concentrations are depleted. We note at least one report of significant difficulty in solving equations with this method, which led the authors to alter the equations at the stage of numerical implementation. 62 For reasons of clarity and simplicity, we recommend that intentionally simple, empirical models should be guided by the Springer or (corrected) Wang-Um models: use a Nernst or linearized Butler-Volmer equation at the anode, and an irreversible Tafel equation at the cathode. These approximations minimize the required parameterization of the kinetic processes, increasing the robustness of computation and facilitating experimental validation.
Wherever simplifying approximations do not apply and a more exact kinetic treatment is required, a rigorous multistep kinetic model should be considered, since the Butler-Volmer equation is not expected to be any more accurate than the simpler models. We consider the recent upsurge in research activity related to the incorporation of such mechanistic models into full PEMFC models to be a positive development; in recommending much simpler kinetic equations as a substitute for the Butler-Volmer equation, we are focused on remedying general trends among modelers focused on practical and consciously lowerfidelity models, and certainly do not wish to challenge progress toward deeper understanding of the catalytic mechanisms. The experimental corroboration of electrode kinetic models and parameterization remains a key challenge, whether the models are chiefly empirical or mechanistically detailed; in particular, we seek to caution against uncritical re-use of equations and their corresponding parameterization between distinct experimental conditions without suitable validation.

Other Comments on PEMFC Kinetic Literature
Our general review of the literature has revealed some other common misconceptions or unusual definitions in the kinetic equations, which we feel are instructive and are worth highlighting here: • Uncritical use of the law of mass action to infer reaction order from molecular stoichiometry -without any consideration of the plausible mechanism or rate-limiting step of a multistep processyields extreme concentration dependences, especially for the proton species. 37,38 • Some authors define a single set of oxidative and reductive transfer coefficients which are then used in the exponential terms for both the anode HOR and the cathode ORR kinetic equations, as though the two chemically distinct reactions were correlated with each other. [63][64][65][66][67][68]

Conclusions
Although the Butler-Volmer equation, as applied in the PEMFC literature, purports to offer a physically meaningful and general description of electrochemical kinetics, it does not apply generally to multistep reactions, such as the electrode reactions (HOR and ORR) in a PEMFC. Because the kinetics of the HOR at the anode are relatively fast, while those of the ORR at the cathode are highly irreversible, we recommend that it is preferable for clarity and simplicity of models to invoke simplifying approximations explicitly, wherever they are judged suitable. For empirical models which avoid mechanistic de-tail of the electrode processes, use the Nernst equation or linearized Butler-Volmer equation (corrected as described above) for the HOR at the anode, and use the Tafel equation for the ORR at the cathode.
Through a review of the literature, we have shown that, notwithstanding the general criticism of the Butler-Volmer equation above, the equation most commonly called the "Butler-Volmer equation" in the PEMFC literature is not the Butler-Volmer equation. Instead, it substitutes concentration-independent polarization in place of concentration-dependent overpotential, yielding an equation whose concentration-dependence cannot be rationalized based on the straightforward interplay of forward and reverse processes in a reversible electrochemical reaction. This mistake, which apparently originates in works from 1991-1992, 35,36 has been repeated uncritically in over 100 subsequent publications and is predominant among contemporary PEMFC models. The same inappropriate substitution also occurs in linearized Butler-Volmer equations, 33 which have been used for the HOR specifically in over 40 publications.
Under non-standard regimes of operation where simple, empirical models cannot necessarily be used, we have demonstrated the risk of appreciable error in predictions from incorrectly formulated electrode kinetics, especially when spatially resolved data are needed. In such circumstances, there is a need for detailed kinetic analysis employing a specific mechanistic model that is experimentally validated under the operating conditions in question: we urge the wider community of full cell PEMFC modelers to heed recent progress in this area. 16,17,44,47 List of Symbols