Interface model of PEM fuel cell membrane steady-state behavior

Modeling works which simulate the proton-exchange membrane fuel cell with the computational fluid dynamics approach involve the simultaneous solution of multiple, interconnected physics equations for fluid flows, heat transport, electrochemical reactions, and both protonic and electronic conduction. Modeling efforts vary by how they treat the physics within and adjacent to the membrane-electrode assembly (MEA). Certain approaches treat the MEA not as part of the computational domain, but rather an interface connecting the anode and cathode computational domains. These approaches may lack the ability to consistently model catalyst layer losses and MEA ohmic resistance. This work presents an upgraded interface formulation where coupled water, heat, and current transport behaviors of the MEA are modeled analytically. Improving upon the previous work, catalyst layer losses can now be modeled accurately without ad-hoc selection of model kinetic parameters. Key to the formulation is the incorporation of water absorption/desorption resistances. The interface model is developed with the consideration of only thru-plane variation, based upon varied fundamental research into each of the relevant physics. The model is validated against differential cell data with high- and low-humidity reactants. The agreement is very good in each case.


List of symbols
Ohmic (referring to total ohmic resistance) Pt Platinum (referring to platinum material) SAT Saturation (saturated state) T Thermal (referring to heat transport) V Void volume fraction (referring to void space) W Water (water dissolved in ionomer)

Introduction
Polymer electrolyte membrane or proton-exchange membrane (PEM) fuel cells are electrochemical energy conversion devices that produce electrical energy from the chemical energy present in hydrogen fuel. Water and fractional waste heat are the byproducts. As the PEMFC operates at temperatures lower than those of other major technologies, it is a candidate for numerous applications. Interest from the automotive industry in low-humidity operation with thin (~ 30 μm) membranes has been noted. The major focus of fuel cell cost reduction and performance improvement strategies is said to be on the issues of (i) heat and water management and (ii) new material development [1]. The basic components of the PEM fuel cell can be explained by means of a cut-away diagram in Fig. 1. A thin membrane-electrode assembly (MEA) separates the anode and cathode flow regions, sandwiched between the porous diffusion media. The left side is the negative, or anode terminal, and the right side is the positive, or cathode terminal. Electrical connection to an external circuit is made via the electrically conductive current collector plates and diffusion media. Both current collector plates typically have gas flow channels that direct the flow of the hydrogen fuel within the anode side, and the oxygen or air oxidizer on the cathode side. Both gas streams are typically pressurized, humidified, and supplied in carefully metered amounts.
The computational fluid dynamics (CFD) approach is used for PEMFC design and simulation. The computational domain can include the flow channels, diffusion layers, and the three MEA regions shown. The CFD approach has been successfully commercialized; however, the computational costs remain quite high. The PEMFC is a multi-scale problem. Simulation of even greatly simplified PEMFC flowfield designs had been reported in the literature benefitting from advances in parallel computing. Costs were driven by the requirements of meshing/discretizing the flow channels (~ 2 × 10 −3 m thickness) simultaneously with the extremely thin catalyst layers (~ 1-2 × 10 −5 m) and membranes (~ 3-20 × 10 −5 m) that comprise the MEA [2]. Significant validation efforts have followed this approach [3,4].
Interface CFD approaches omit the thin MEA from the computational domain, treating the MEA as an interface, which separates the anode and cathode flow domains [5]. The MEA can be treated as a reacting wall, with consumption/production source terms on either side to mimic its operation. This approach may entail less computational cost. It has been criticized in the literature, however, for lacking the capability to model reaction rates and losses occurring within the thin, porous catalyst layers, detailed water distribution in the membrane, and the important transient effects linked to these phenomena.
Since these interface approaches were first published, much research has improved the understanding of the multiple MEA physics. The objective of this work is to reformulate and improve the interface approach to represent the MEA, accurately modeling the relevant physics (including catalyst layer losses), for steady-state problems. Improvements are meant to (i) better model the underlying MEA physics based upon the best available published results, and (ii) improve the methods for devising these calculations. Flow equations from the 3-D computational domains are not the focus here. Validation against differential cell data is undertaken where gas composition is known/uniform. The MEA composition must be thoroughly described for successful validation.

Model development
The MEA is represented as a two-dimensional interface that separates the anode and cathode computational domains. The MEA computational routine accesses the solution variables at the surface of each domain and produces fluxes into each of the physics. The role of the interface model in multi-domain approaches has been to approximate all the externally relevant behavior of the MEA (i.e., current generation, reactant consumption, water permeation, heat generation, etc.) by boundary conditions on both sides of the interface [6]. Figure 2 shows an MEA schematic with typical thickness dimensions included. The anode contains hydrogen, water vapor, and inert nitrogen species, while the cathode is modeled as a mixture of oxygen, water vapor, and inert nitrogen. The interface takes as boundary values, the anode (A) and cathode (C) solution variables from the GDL. The interface replaces the MEA, containing the anode catalyst layer (ACL), membrane (MEM) and cathode catalyst layer (CCL) The MEA model has inputs and outputs shown in Fig. 3. Inputs include gas pressure, temperature, mole fractions, and voltage from the adjacent anode and cathode cells within the gas diffusion layer (GDL). Mole fractions of the species present, such as water, hydrogen, and oxygen, are used similarly to the previous works. The outputs include estimates of current density, heat sources, reactant fluxes, water fluxes, and the water gain rate of the ionomer phase of the MEA.

MEA composition/description
A description of MEA dimensions and compositions is needed to estimate the various catalyst layer losses. The MEA has three regions: the anode catalyst layer, membrane, and cathode catalyst layer, denoted by ACL, MEM, and CCL, respectively.
The membrane is assumed to be Nafion, a solid ionomer which is nearly impervious to gas penetration, except for water vapor. A constant, selectable equivalent weight EW (g/equiv or g/mol SO 3 − ), and density Io (g/cm 3 ), is contemplated. Absorbed water with density W (g/cm 3 ) exists within the ionomer. A constant crossover current density, I x (A/m 2 ) is assumed: typically, 0.5-1 mA/cm 2 (5-10 A/ m 2 ), as a result of slight hydrogen gas permeation through the membrane from anode to cathode. The thicknesses of the three MEA zones, denoted t ACL , t CCL , and t MEM (m), refer to the dry values; prior to water uptake. Nominal membrane thicknesses refer to values at about 50% RH; hence, a 22 μm-thick dry membrane swells to a 25 μm nominal thickness [7].
PEMFC catalyst layers are porous electrodes: with mixtures of platinum catalyst nanoparticles, the carbon support, the ionomer binder, water sorbed within the ionomer, and void space, that allow for reactant gas diffusion. The ionomer is the same EW as the membrane. The material components, their respective functions within the catalyst layer(s), and their transport roles are summarized in Table 1. To define catalyst layer composition, several densities, molecular weights, and fundamental physical constants are listed in Table 2. Catalyst layer composition described in Table 3 is assumed to be spatially uniform throughout its thickness. The composition of the anode catalyst layer (ACL) and cathode catalyst layer (CCL) is determined by four parameters: the platinum loading L ACL,Pt (~ 0.1-0.4 mg Pt cm −2 ), the weight percent Pt/carbon in the catalyst powder Ptc ACL (~ 20-60%), the ionomer-to-carbon ratio IC ACL (~ 0.5-1.5) which is the mass ratio of dry ionomer to carbon within the respective catalyst layer, and its thickness t ACL (~ 2-12 μm). Once these are known, the composition of the catalyst layer is further described.

Ionomer water uptake
Ionomer equilibrium water uptake is expressed in the form of a non-dimensional water uptake eq (mol H 2 O/mol SO 3 − ) the number of water molecules per acid site within the   ionomer, which normalizes water uptake for any equivalent weight. Uptake from humidity and ionomer swelling with water uptake is described. Ionomer water uptake is modeled with reference to water activity, or gas humidity, at anode ( A ) and cathode ( C ) side adjacent cells using ideal gas mixture relationships. Relative humidity is calculated as in Eq. (1), where water vapor saturation pressure P SAT (T) is defined by a fit in Eq. (2) with temperature (K) recommended by [8]. Equilibrium water uptake is curve-fitted from the experimental measurements of Jalani et al. [9,10] with fifth-order polynomials. Uptake from Nafion membranes in contact with liquid water has been observed at eq = 22 water molecules per acid site, while uptake from the vapor phase ranges from eq = 1 − 2 (10% RH) to eq = 14 − 16 (100% RH). Water uptake from the liquid phase (300% RH) is eq = 22 , in agreement with several previous authors [11,12]. Detailed derivation can be found in Edwards [13]: Water uptake within the ionomer will vary from the equilibrium values eq . The ionomer is considered a mixture of the two phases: dry membrane and water. The concentration of dissolved water in the ionomer, c W (mol H 2 O m −3 ), is related to the dimensionless water content expressed as in Eq. (3), where c f is the concentration of fixed charge sites (mol SO 3 − m −3 ). The volume fraction of water within the ionomer phase, f V , is given as Eq. (4), where W represents the molar volume of water (m 3 mol −1 ) and Io the Ionomer volumetric swelling results in increased MEA thickness. While the in situ PEMFC environment compresses the MEA, a focused study concluded that, at typical cell assembly pressures, the water uptake is not significantly decreased [14]. The swollen membrane thickness with water uptake is given in Eq. (6), similar to the experimental kinetics investigations of Liu et al. [7]. Though not as obvious, those investigations contemplate a similar expansion in catalyst layer thickness:

Ionomer-phase water/ion transport
The membrane is modeled as impermeable to gases, but considers the absorption, desorption, and permeation of water.
The current density ⃗ Io (A m −2 ) represents protons (H+) moving through the membrane under the influence of the gradient of the electric potential field of the ionic phase Io . Consumption of hydrogen and oxygen reactants is accounted for by boundary fluxes. Water and current transport within the membrane are assumed to occur only in the thru-plane direction and are treated with the common, de-facto standard phenomenological approach originated by Springer et al. [15], using only macroscopic calculations. The flux forms are shown in Eq. (7). When restricted to the thru-plane direction, current density ⃗ Io = I must be constant through the thickness of the membrane: The three transport properties n d , D W , and Io exhibit water content and temperature dependence. Property variation makes this problem inherently non-linear and the degree of property variation with changes in water content was said to be substantial [16].
The electroosmotic drag coefficient, n d (-), represents the number of water molecules dragged from anode to cathode along with each migrating proton. Typical approximate values have been reported as 0.6-1.0 for vapor-equilibrated membranes and 2-2.5 for liquid-equilibrated membranes [17]. The diffusion coefficient of water in the Nafion membrane phase, D W (m 2 s −1 ), is also a complex function of temperature and water content. This work will use the methodology of Ge et al. [18], where the diffusion coefficient was found to be between 2 and 10 × 10 −10 m 2 s −1 , varying with water content and temperature. The ionic conductivity of the ionomer Io , (S m −1 or Ω −1 m −1 ), is also very much water content and temperature-dependent. Typical values range from 2 to 20 S m −1 , with the greatest conductivity being in liquid-equilibrated membranes [19][20][21].
Interfacial resistance to water transport must be considered. Water absorption into the ionomer, and desorption out of the ionomer phase occurs as Eq.
a net normal water flux magnitude. The terms k a and k d represent mass-transfer coefficients for the absorption and desorption, respectively, of water to and from the ionomer. They are themselves dependent upon the membrane water content: Relations for the transport properties n d , D W , Io , and k a |k d are needed. The relations follow the approach of Ge et al. [17,18]. However, the 1994 water uptake curves of Hinatsu et al. [22], used in that work, are replaced with the more recent 2005 water uptake curves of Jalani and Datta = eq (liquid water).
[10]. Updated conductivity relations are included, as well. Further details are available in Edwards [13]. These newer water uptake data are believed to be more representative of the in situ PEMFC environment. The ionic conductivity of Nafion 11 × membrane at 80 °C and 100% humidity is known to be around Io = 17 Ω −1 m −1 [19]. The water uptake curve of Hinatsu et al. [22] predicts a water uptake of = 9 under these conditions; when combined with the popular conductivity relations of Springer et al. [15], it gives Io of only 8 Ω −1 m −1 , creating an obvious discrepancy. Later work by Onishi et al. [23] in 2007 explicitly claimed that the 1994 water uptake data were biased. Examining the presence of "Schroeder's paradox" in Nafion water uptake measurements, they suggested that a pre-drying procedure, at high temperature under vacuum, may have inappropriately influenced the membrane's morphology, and reduced the water sorption of the 1994 test. Therefore, the transport relations used in the prior work of Ge et al. [17,18] are recomputed incorporating the 2005 water uptake data of Jalani and Datta [10].

Ionomer-phase water content profile
An MEA hydration and temperature profile are needed. Approximate heat and water transport equations must be solved within the MEA to estimate ohmic resistances, effective catalyst layer losses, and current density. The interface model iteratively calculates (i) a water content profile, (ii) current density, and (iii) a temperature profile. Figure 4 shows a profile schematic to explain the water and heat transport equations. An approximate water content profile (z) is assumed to be a second-order polynomial determined by three unknown values ( 1 2 3 ) occurring at the ACL-membrane interface (anode side), midpoint, and CCL-membrane interface (cathode side) of the membrane. The equilibrium water uptake of the anode, eq A , and cathode, eq C , are found using the water activity and temperature of the anode and cathode gases. The water contents of the ionomer in the ACL, ACL and CCL, CCL , are assumed uniform throughout the respective catalyst layers, with the values of The water transport model employs convective boundary conditions of water absorption and desorption. The MEA can gain or lose water content. The water content profile utilizes three equations: two values of its slope, at points 1 and 3, in addition to the integrated water content W , at pseudosteady-state conditions. The pseudo-steady-state assumption involves ignoring water-content distribution variations within the ionomer of the MEA, to represent hydration level by an overall MEA water content.
The water contents are found by assembling three equations to solve for the three unknown water contents. The first equation relates the MEA water content W to the three water content values. The MEA water content W (mol m −2 ) is represented as a lumped parameter value for each cell of the interface. During steady-state operation, the value of W must be such that the water gain rate W∕ t = 0 . When W is less than the steady-state value, the water gain rate is positive; the converse occurs, as well. The overall water content of the ionomer phase, integrated through the thickness of the MEA can be presented as Eq. (9), which simplifies to Eq. (10): The second equation examines the slope in the water content profile at z = 0, the ACL-membrane interface. There is water flux to/from the ionomer by absorption or desorption from the anode gas stream, according to the absorption/desorption kinetics previously given. The water flux within the membrane should be identical, being driven by the combination of electroosmotic drag and back diffusion. The absorption flux (positive) or desorption flux (negative) is presented as Eq. (11) and the water flux within the ionomer phase at z = 0 is written as Eq. (12) using a three-point estimate of the slope at z = 0: The expression simplifies to Eq. (13). The transport properties (drag coefficient n d,1 , diffusion coefficient D W,1 , and mass-transfer coefficient k a|d ) are evaluated at the water content and temperature of the ACL-membrane interface The third equation sets the slope in the water content profile at z = t MEM , the CCL-membrane interface. There is a water flux within the membrane driven by electroosmotic drag and back diffusion. Product water is produced in the ionomer phase of the CCL by the ORR. To the right of the CCL-membrane interface, there is desorption or absorption flux into or from the cathode gas stream. The desorption flux (positive) or absorption flux (negative) is presented as Eq. (14) where the water volume fraction term uses 3 as the water content. The water flux within the ionomer phase at z = t MEM is written as Eq. (15) using a three-point estimate of the slope at z = t MEM : Water production by the ORR also occurs here. If interfacial resistance to water transport is incorporated into modeling efforts, a decision needs to be made as to where water creation occurs: in the CCL gaseous phase, ionomer phase, or directly as liquid. The extensive investigations of Wu et al. [24,25]  have explained that water creation takes place in the ionomer phase.
The flux into the cathode, by electroosmotic drag and diffusion, is added to the product water creation, and set equal to desorption/absorption flux exiting the cathode into the cathode gas stream. Equating the water fluxes J W,C yields Eq. (16), where k 3 is evaluated as an absorption or desorption coefficient according to the relative values of eq C and 3 . This expression simplifies to Eq. (17) where, again, the coefficients of diffusion ( D W,3 ) and electroosmotic drag ( n d, 3 ) are evaluated at the water content and temperature of the CCL-membrane interface ( 3 , T 3 ): Together, there is a system of three equations and three unknowns to solve for the water contents in Eq. (18). The solution updates ( 1 2 3 ) . The current density ( I ) is updated in a subsequent step. After the water content profile is calculated, the fluxes of water into the anode, J W,A , and out of cathode, J W,C , (mol m −2 s −1 ) can be updated from the water content profile of the MEA as Eq. (19): The water gain rate within the MEA, W∕ t (mol m −2 s −1 ), is found by applying conservation of mass principles to the water in the ionomer phase. In Eq. (20), the first term represents water production due to the ORR; which takes place in the ionomer phase. The next term represents water flux into the anode and the final term water flux out of the cathode:

Gas-phase transport in catalyst layers
Gas transport losses within the anode and cathode catalyst layers are considered in this section. Hydrogen and oxygen (16) reactants are consumed through the thicknesses of anode and cathode porous catalyst layers (electrodes). Reaction rates will vary with depth through the catalyst layer; however, the distribution of reaction rates is predominantly shaped by conduction losses in the catalyst layer's ionomer [26,27] and not by normally occurring variations in reactant gas concentration. The reactant mole fractions at the outer edge of the ACL and CCL, X H 2 ,A , and X O 2 ,C are adjusted to diffusioncorrected values X H 2 ,A and X O 2, C at the location of average reaction current density as shown in Eqs. (21) and (22). The normal current density is I (A m −2 ); t ACL and t CCL (m) represent the effective gas diffusion lengths within each catalyst layer, and D eff H 2 ,ACL , D eff O 2 ,CCL are the effective catalyst layer diffusion coefficients: The through-plane gas-phase transport of the reactants occurs due to diffusion; pressure-driven flow is not considered. The consumption of reactants causes a gradient in the reactant mole fractions between the outer edge of the ACL/ CCL and the membrane-catalyst layer interfaces. A simplified approach is employed; it assumes that all the reactants to be consumed (by stoichiometric requirements) must diffuse over an effective length. Effective length is treated as a known here, but it is determined analytically by conduction effects in a later section. The calculation details of the effective catalyst layer diffusion coefficients were given in [13]. Both molecular diffusion and Knudsen diffusion mechanisms must be accounted for. Porosity and tortuosity effects can be accounted for to produce an effective oxygen or hydrogen diffusivity. A reduction in mole fraction, at the catalyst layer's mean reaction depth, can be calculated given the reactant flux necessary for the given current density.

Kinetics in catalyst layers
Kinetic losses, or activation losses, within the anode and cathode catalyst layers are considered in this section. These are the voltage losses required to drive the chemical reactions at the rate (current density) required. Anode kinetics typically is not that important to a PEMFC model as the hydrogen oxidation reaction (HOR) is quite facile, with a large exchange current density. Cathode kinetics is very important to PEMFC modeling and the subject of a great deal more research. The oxygen reduction reaction (ORR) occurring at the cathode is a source of major voltage loss with currently available practical PEMFC catalyst configurations. The kinetic models and parameters associated with each reaction need to be detailed. These are gathered from focused research in each area. This information is then subsequently used to calculate the ohmic loss terms within the MEA and estimate current density.

Anode kinetics
Anode kinetics follows the general Butler-Volmer (BV) model. The ideal BV model pre-supposes an anode catalyst layer (ACL) without ohmic resistance or gas diffusion losses [8,28] in Eq. (23), where I (A m −2 ) refers to the cell current and I 0,A refers to the superficial exchange current density for the HOR occurring at the anode electrode (explained below). Here, the term A refers to the kinetic or activation overpotential at the anode (V): the voltage loss. The terms a and c are the anodic and cathodic charge transfer coefficients, respectively. These are fundamental constants describing the reaction, which are not supposed to change with operating conditions such as temperature, humidity, etc. [29,30]. They are not used here as adjustable parameters in a given problem. The superficial exchange current density I 0,A is the dominant parameter in the BV equation to be influenced/optimized, influencing the voltage loss. Exchange current density depends on gas temperature, pressure, and reactant concentration (mole fraction). It also depends on electrode morphology, the catalyst type, and the aging/degradation of the catalyst layer. It is known to be an exponentially increasing function of temperature: The HOR reaction is facile; the exchange current density can reach ~ 20 times or more the size of the maximum operating current density at typical operating conditions. The bracketed term must necessarily be small as well and so the arguments to the exponential terms are typically less than 0.15. With this limitation, the linearized BV equation in Eq. (24) can represent the losses of the HOR [8,28]: Several sources indicate that the reaction is symmetric, i.e., a = c [28,[30][31][32]. Sources indicate that a = c = 0.5 , and thus, a + c = 1 [8,28,30]. Other sources indicate that a + c = 1 , and thus, a + c = 2 [31,32]. Focused experiments [29] claimed that both values of the charge transfer coefficients could fit the experimental data. As of 2015, Durst et al. showed that a + c = 1 could fit measurements performed over a range of temperatures [30].
The superficial exchange current density term, I 0,A (A m −2 ), is a product of several factors in Eq. (25); the form attempts to separate out and identify the different influences of temperature, gas concentration, catalyst layer morphology, and platinum loading. The exchange current density is proportional to the mole fraction of hydrogen X H 2, A present, at the mean reaction depth, in the anode catalyst layer. The term i 0,ACL ( A cm −2 Pt ) refers to the specific exchange current density of the HOR on platinum catalyst, expressed with the basis of the active area of platinum catalyst. It was reported [30] to be between 0.24 and 0.26 ( A cm −2 Pt ) at a temperature of 80 °C, and a temperature dependence was investigated. L ACL,Pt (mg Pt cm −2 ) is the anode platinum catalyst loading, known from the catalyst layer composition: The catalyst-specific exchange current density for the HOR on the carbon-supported platinum catalyst, i 0,ACL ( A cm −2 Pt ), has dependencies of temperature and hydrogen partial pressure in Eq. (26). The term i * 0,ACL represents the value at reference hydrogen partial pressure of 101 kPa and 80 °C (353 K) of 0.26 A cm −2 Pt [29,30]. The activation energy of the HOR, E rev A (kJ mol −1 ), was found to be 16 kJ mol −1 . The dependence of anode losses on hydrogen concentration is modeled through the diffusion-corrected hydrogen mole fraction X H 2, A present at the mean reaction depth in the ACL. Mass transport or diffusional losses result when H 2 molecules cannot reach the CCL reaction sites to support the desired reaction rate. As X H 2, A diminishes, superficial exchange current density I 0,A decreases and the anode overpotential A grows:  The term A ACL,Pt ( m 2 Pt ∕g Pt ) is the anode catalyst layer's electrochemically available catalyst area, (ESA or ECSA), expressed per gram of catalyst, indicating morphology effects. This term measures the effectiveness of the catalyst layer microstructure and design in employing the catalyst to create maximum effective active area where the reaction can take place. It represents the surface area of the Pt catalyst that is simultaneously in contact with the reactant gas, the electronic conducting phase, and the proton conducting phase of the catalyst layer. It has been found to be typically 46-60 m 2 Pt ∕g Pt [7,8] at fully humidified conditions. Experiments with typical catalyst layer structures showed that this area could drop by 50% as humidity decreases from 100% RH to 50% RH [7,[33][34][35].
Carbon black supports, such as the commercially available product Vulcan XC-72 (Cabot Corp), were seen with relatively large pore-size distributions (PSD). Carbon supports with larger PSD tend to have higher ionomer coverage [34] of the catalyst sites, and less reliance on surface-adsorbed water for proton conduction within the CL. For such a carbon black support, the ECSA is treated here as decreasing linearly with local CL water content , and not humidity [7], from 0 to 100% RH in Eq. (27). This introduces more nonlinearity into the problem of determining the current density and water contents within the MEA: This section described the HOR kinetic losses and superficial exchange current density I 0,A (A m −2 ), occurring at the anode electrode. The purely kinetic losses of the anode are typically small and have been routinely ignored in modeling efforts [36]. Ohmic losses in the ACL have not yet been considered; these can become significant under sub-saturated conditions. Combined ohmic and kinetic losses are needed, and will be addressed subsequently.

Cathode kinetics
BV kinetics simplify to Tafel kinetics when modeling the ORR occurring at the CCL. With very low exchange current density, and without significant conduction or diffusional losses, Tafel kinetics follows Eq. (28) [37] where the term C (V) refers to the purely kinetic or activation overpotential of the ORR at the cathode. I (A m −2 ) refers to the cell current as before and I x the crossover current density (which represents a small amount of hydrogen permeation "crossover" from the anode through the membrane into the cathode). The term I 0,C refers to the superficial exchange current density for the ORR occurring at the cathode electrode. Exchange current density for the ORR taking place on platinum catalyst depends on gas temperature, pressure, reactant concentration (mole fraction), electrode morphology, and the catalyst type. It is known to be an exponentially increasing function of temperature. The term c is the cathodic charge transfer coefficient, which was found experimentally to be unity for the ORR when ohmic losses were correctly accounted for [37]. The Tafel slope indicated is 70 mV/decade at 80 °C (with c = 1): Experimental works have shown considerable disagreement in measuring the Tafel slope. Some works have reported "double Tafel slopes" where the values consistent with c = 1 are reported at low-current densities together with a doubled Tafel slope ( c = 0.5 ) at greater current densities. Prior CFD models, utilizing the interface approach, have employed doubled, or nearly doubled, Tafel slopes [38,39] for use with low-humidity conditions. In the present model, the cathodic charge transfer coefficient c = 1 and does not change with current level, humidity, or other factors.
The superficial exchange current density for the ORR I 0,C (A/m 2 ), occurring at the cathode electrode, is a product of several factors in Eq. (29); this form attempts to separate the different influences of temperature, gas concentration, catalyst layer morphology, and platinum loading. I 0,C is proportional to the mole fraction of oxygen X O 2, C present, at the mean reaction depth, in the anode catalyst layer. The term i 0,CCL ( A cm −2 Pt ) refers to the catalyst-specific exchange current density of the ORR on platinum catalyst, expressed with the basis of the active area of platinum catalyst. The term L CCL,Pt ( mg Pt ∕cm 2 ) is the cathode platinum catalyst loading, known from the catalyst layer composition, and A CCL,Pt ( m 2 Pt ∕g Pt ) is the cathode catalyst layer's electrochemically available surface area (ECSA), which is as previously discussed for the anode: The catalyst-specific exchange current density for the ORR on the carbon-supported platinum catalyst, i 0,CCL ( A cm −2 Pt ), has dependencies of temperature and oxygen partial pressure in Eq. (30) where i * 0,CCL represents the value at reference oxygen partial pressure of 101.3 kPa and 80 °C (353 K). It was found to be 2.5 × 10 −8 [37] and later 1.5 × 10 −8 A cm −2 Pt [7]. The kinetic reaction order with respect to oxygen partial pressure was found to be 0.54, and E rev C , the activation energy of the ORR (at zero overpotential/ reversible cell potential), was found to be 67 kJ mol −1 :

Ohmic losses within the MEA
This section details the effective ohmic resistances present within the MEA. These resistances are needed for the current estimation. Ohmic losses within the interface are calculated from the hydration-dependent conductivities.
In situ high-frequency resistance (HFR) measurement produces R (Ω m 2 ) that is the sum of the membrane protonic resistance and the electronic resistances present as in Eq. (31). The HFR measurement is commonly used for the model validation of water transport modeling efforts, because the HFR measurement contains a first component that has membrane water-content dependence and a second component that has been found to be constant for a given cell assembly. The total ohmic resistance within the interface, R Ohm (Ω m 2 ), is the sum of HFR resistance and effective catalyst layer resistances as Eq. (32):

Effective electronic resistance
The Interface model MEA contains electronic conduction resistances R e − cnt (Ω m 2 ). This term is treated as a humidityindependent, but assembly-pressure-dependent, term. It can represent ACL/CCL electronic conduction losses and contact resistances between the gas diffusion media (GDL) and ACL/CCL (catalyst layers) at the outer edges of the MEA.
When comparing the interface model directly to experimental data, additional bulk resistances must be included in the electronic resistance term. Electronic bulk conduction losses may occur outside of the interface (MEA) region; in the anode and cathode diffusion media and the conductive flow-field plates. Electronic contact resistances can also be found between the diffusion media and flow-field plates.

Membrane resistance
The membrane proton conduction resistance R H + MEM (Ω m 2 ) is calculated by an approximate integration of the three hydration-dependent conductivity values. Volumetric swelling in the membrane is treated as occurring purely in the thruplane direction (in thickness); membrane thickness swells with water content at the center at the membrane as following equation:

Effective anode resistance
The voltage drop occurring in the ACL is modeled by an analytical solution which accounts for both kinetic and ohmic losses simultaneously. Kulikovsky [28] developed an analytical solution for anode overpotential, in the absence of gas-phase diffusional losses (assuming ideal reactant transport), which will be expanded somewhat to account for ACL swelling which was not originally considered.
The ACL overpotential first depends on the ACL electrode sheet resistance, R H + ,sheet ACL (Ω m 2 ). The sheet resistance of the anode is the protonic conduction resistance through the electrode, assuming that all current travels through the entire thickness of the electrode. The sheet resistance calculation in Eq. (34) assumes that the ACL swells in thickness with water uptake. The term t ACL (m) represents the CCL dry thickness, ACL,Io (-) the dry volume fraction of ionomer in the CCL, dV,1 (-) the fractional volume swelling with water uptake at the anode side of the MEA, Io,1 (Ω −1 m −1 ) the ionic conductivity calculated from water content and temperature, and ACL, Io (-) the tortuosity of the ionomer within the ACL. The tortuosity of the ACL ionomer conduction network is calculated in the same manner as for the CCL: The anode catalyst layer thickness can be viewed as being split into a conduction and diffusion thickness. The effective conduction thickness, t ℂ ACL , and the effective gas diffusion thickness, t ACL , sum to the swollen anode catalyst layer (ACL) thickness in Eq. (35). The HOR reaction occurs within a thin strip at the membrane side of the ACL, and so, the ACL diffusional thickness is essentially the entire catalyst layer thickness. This effective diffusion length is used in Sect. 2.6 to ascertain the reduction in reactant mole fractions within the anode catalyst layer: The anode loss A (V) and the operating current density I are related as Eq. (36), where eff ACL is the catalyst layer effective conductivity. The quantity under the radical is made small (< 1) by the large superficial exchange current density of the HOR, I 0,A , making the sinh function nearly linear in character with respect to current density I . To a first approximation, overpotential A can vary linearly with current. The term I 0,A will drop with diffusional reductions within the ACL, and with reductions in the ACL ionomer-phase water content (due to electroosmotic drag). The ACL conductivity eff ACL will also decline with increasing current density for the same reason. Accounting for the thickness expansion of the ACL, this formula can incorporate the effective anode sheet resistance R H + ,sheet ACL of the ACL as Eq. (37): If a small argument approximation is applied to linearize the arcsinh function, the relation can be re-written to appear explicitly linear as Eq. (38). The ACL voltage losses have been assessed, experimentally, as being not completely linear with current density in a hydrogen concentration cell [19,29,30]. This last linearization is not employed here. The effective catalyst layer conduction resistance R H + ,eff ACL (Ω m 2 ) is calculated from Eq. (39), where the combined kinetic and ohmic losses of the ACL in Eq. (38) are linearized to an effective resistance: Finally, the bounds of validity of this anode loss approximation are addressed. This anode loss relationship was formed using a small argument simplification, described as the low-current approximation. The anode loss estimate has Eq. (40) as a restriction on its validity: With practical ACL compositions, the right-most ratio is significantly greater than 1, and so, the coth function produces a result of very nearly unity. This relationship is thus valid when the ratio of ACL kinetic resistance R H + ,kinetic ACL to ACL sheet resistance R H + ,sheet ACL is suitably small in relation to the current density. It was noted that the PEMFC anode consistently obeys this low-current approximation, termed as the linear mode [28]. Should ACL catalyst loading be greatly reduced, the validity of this relation needs to be re-assessed in terms of the reduced superficial exchange current density I 0,A to be expected [40].

Effective cathode resistance
The effective cathode catalyst layer resistance is derived from an analytical solution that was presented by Wang and Feng [26,27], as well as Neyerlin et al. [19,37].
The analytical solution was developed incorporating several simplifying assumptions within the catalyst layers. The spatial distribution of all material phases (volume fractions) within the CL is assumed constant, as is the ionomer water content within the CCL, where there is constant ionic conductivity. Another key assumption of this analytical solution is that the solid phase, or electronic phase, is thought to have constant phase potential through the thickness of the electrode; electrical conduction losses are assumed insignificant. Electrical contact resistances were incorporated at the outer edge of the catalyst layers. The CCL is assumed to be effectively isothermal [26]. Under these conditions, variations in reaction rate within a porous electrode are dependent upon conduction losses and the Tafel slope of the reaction.
Reactant gaseous concentration levels are assumed to be constant within the variable depth of the catalyst layer, which is also referred to as the absence of O 2 diffusion resistances. Wang and Feng estimated that, in most ranges of fuel cell operation, diffusional losses within the CCL would be insignificant [26]. Detailed numerical analysis of the reaction rate distribution throughout the thickness of a PEMFC catalyst layer has also been performed by modeling. Wu et al. carried out the full analysis numerically, and found that reaction rates located to minimize conduction losses, without noticeable impact from gas concentration losses [41]. This was said to be due to the form of the Butler-Vollmer equation, where the variations in the exponential term (for overpotential) had much more impact on reaction rates than the variations in the linear term (where exchange current density carries a linear dependence on the gaseous reactant concentration term). This work utilizes this analytical solution, but with the small diffusional correction of Sect. 2.6 applied.
The electrode sheet resistance R H + ,sheet CCL (Ω m 2 ) in Eq. (41) expresses the effective protonic conduction resistance through the electrode, assuming that all current travels through the entire thickness of the electrode, where the electrode swells in thickness with water uptake as the membrane does. The term t CCL (m) represents the CCL dry thickness, 1 3 and CCL,Io (-) the tortuosity of the ionomer within the CCL. Details about the calculation of tortuosity are given in [13]: A dimensionless performance parameter CCL in Eq. (42) is the ratio of proton conduction resistance to kinetic resistance within the CCL [19]. The numerator includes the sheet resistance and current density. The denominator is the Tafel slope of the ORR expressed as b in V/decade. It has been found as 70 mV/decade at 80 °C [37]: The effective catalyst layer conduction resistance R H + ,eff CCL (Ω m 2 ) is then calculated from the sheet resistance R H + ,sheet CCL and the resistance correction factor as CCL ( CCL ) as Eq. (43) [19]. The correction factor CCL ( CCL ) is fit with an approximate function in Eq. (44) valid between 0 and 30. Values below and above that range were limited to 1 and 5.908, respectively: The thickness of cathode catalyst layer can then be viewed as being split into a conduction and diffusion thickness. The effective conduction length, t ℂ CCL , and the effective gas diffusion length, t CCL , are defined by Eqs. (45) and (46). This effective diffusion length is used in Sect. 2.6 to ascertain the reduction in reactant mole fractions within the cathode catalyst layer:

Polarization curve and current estimation
This section describes the voltage-current relationships in the interface model. Established voltage and current relationships are used to produce an updated estimate of current density with the Newton-Raphson technique.
The thermodynamically determined open circuit, ideal, or reversible voltage can be calculated from the established thermodynamic relationships. The relationship used in this work was given by Liu et al. [7]. The reversible cell voltage  (45) t ℂ CCL = t CCL 1 + dV,3 .
V OC is Eq. (47). The reversible cell voltage can be seen to be a function of cathode temperature ( T CCL ), hydrogen partial pressure ( X H 2 ,A P A ), oxygen partial pressure ( X O 2 ,C P C ), and cathode humidity C . The definition of the open-circuit voltage takes into account some reduction in the oxygen and hydrogen mole fractions due to mass-transfer (diffusional) losses within the catalyst layers: The cell voltage in Eq. (48) is the difference between cathode ( V C ) and anode ( V A ) potential levels. Cell current is the previously described current density. The cell voltage is the ideal/reversible voltage minus various voltage losses [19]. The ACL combined kinetic and ohmic losses have been incorporated into R H + ,eff ACL within the ohmic resistance term R Ohm . The iR-free or kinetic voltage V k in Eq. (49) describes the cell voltage in the absence of ohmic losses. The kinetic voltage is formed analytically by subtracting cathode kinetic losses C from V OC : The voltage loss V loss is defined as a simple re-arrangement to give Eq. (50), which is the voltage-current relationship of the interface, where the current density I needed to satisfy this relationship is sought, but is not known. An implementation of the Newton-Raphson root-finding method to solve for the current density I new was found sufficiently stable for use here. The procedure initializes I = I x and then calculates an updated estimate of the operational current density I new using Eq. (51):

MEA thermal treatment
The interface takes as inputs the temperatures from the opposite sides of the MEA and outputs heat fluxes to the anode and cathode computational domains. Some earlier interface models created a thermal source term at the cathode GDL/MEA interface based upon the (estimated) overall efficiency of the device [42]. This present model divides the heat flux between anode (47) and cathode, but does not go into high detail examining the spatial distribution of heat generation within catalyst layers. It is included, because several terms from the interface model were given temperature-dependent. A 1-D temperature profile estimates the temperature rise in the MEA in Fig. 4. Temperatures T A and T C are taken as Dirichlet boundary conditions from the 3D simulation as thermal contact resistance is thought to be minimal. A simultaneous solution of the seven temperature values (five unknowns) is developed to produce a steady-state temperature profile and heat fluxes into the anode and cathode flow domains.
This section applies a general heat transport equation to the MEA. Heat generation within the regions of the MEA has been detailed [43]; it typically occurs predominantly within the CCL. The thermal conductivities and heat generation terms of each region of the MEA are collected. Convective effects within the MEA are ignored, following similar work [44]. Heat transport in the thru-plane direction occurs by conduction, with the thermal conductivity being hydration-dependent. The regions of the MEA are thought to be isotropic and thermal contact resistances between the layers of the MEA are thought to be negligible [45].
The governing equation of heat transport can be written as Eq. (52), where eff i refers to the thermal conductivity (W m −1 K −1 ), T the temperature (K), Q i the superficial heat source (W m −2 ), and t i (m) the layer thickness. The subscript i refers to the three regions of the MEA (i = ACL, MEM, CCL). Details of thermal conductivities, eff i , of the MEA regions are given in [13]: Estimates of the superficial heat generation, Q i , within the regions of the MEA, are available [24,44]. Within the anode catalyst layer, heat generation occurs from irreversible heating sources. The combined kinetic and ohmic ACL losses will generate heat, in addition to half of the MEA electronic resistance. The anode superficial heat generation is expressed as Eq. (53) where the ACL swells in thickness with water uptake: Within the membrane, there is only the irreversible heat source of ohmic losses (protonic conduction losses). The membrane superficial heat generation is expressed by the following equation: The cathode catalyst layer has both reversible and irreversible heat sources in Eq. (55). Reversible heat generation (53) .
in the CCL stems from the entropy change of the overall reaction, and irreversible heat generation from the cathode activation (kinetic) and ohmic losses. The entropy change of the reaction, dS (= 44 J mol −1 K −1 ), assumes that produced water exits the MEA in vapor form [8,24,32], and is approximately constant at temperatures below 100 °C. The term C represents the purely kinetic losses in the cathode as described earlier. The right-most term represents ohmic heat generation, utilizing the effective cathode resistance and half of the MEA electronic conduction resistance: The temperature profile within the MEA is now to be determined. Referring to Fig. 4, there are seven temperature values in the approximate profile, and thermal contact resistances are negligible [45]. The regional heat sources are approximately uniform (within each region). The temperature profile can be expressed as the solution to a system of five heat transfer equations with the five unknown temperatures.
The conservation of energy equation, Eq. (52), is applied to the anode catalyst layer in Eq. (56), the membrane in Eq. (57), and the cathode catalyst layer in Eq. (58): heat fluxes into the anode J T,A and cathode J T,C domains can be solved for in Eq. (63), again using three-point approximations to the first derivative of the temperature profile:

Boundary conditions
When the iterative scheme of the previous sections converges on a MEA current density estimate, temperature profile, and water-content profile, the interface model produces its outputs to the 3-D computational domains. Boundary conditions on the anode and cathode side of the MEA represent current density, heat flux, reactant consumption, and water flux in Eq. (64). Water content is also updated:

Iterative operation
This section describes the iterative operation of the interface model. As described, this algorithm can be run within a general-purpose computing package such as MATLAB. The main variables of the iterative routine are given in Table 4. In the inner loop, the water content W and water fluxes J W,− are held constant, when the remaining variables are solved for.
First, the water content profile 1 2 3 is determined. Then, volume fractions within the ionomer are calculated. Water transport parameters are then determined for the ionomer phase. The water transport equations are then solved with the most recent estimates of current density J e,A = I J e,C = −I and the other solution variables. Second, the MEA ohmic resistances are re-assessed. Effective membrane, anode, and cathode resistances are updated from the water content profile. Third, gas-phase diffusional adjustments, representing mass-transfer resistances within the catalyst layers, are found. The diffusion-corrected oxygen and hydrogen mole fractions are then calculated. Fourth, the current density estimate I is updated. Finally, the MEA temperature profile is updated. The outer loop operates to adjust overall hydration and boundary water flux values, in addition to the remaining boundary conditions described in the previous section. These calculations run only after a converged estimate of current density I is achieved within the inner loop. Water fluxes J W,A and J W,C are updated, but with relaxation factor (~ 0.5) to ensure stability and convergence. Water content W is updated with a pseudo time step (~ 0.25-0.5 s). The water gain rate W∕ t will approach zero as W nears its equilibrium value. The convergence of the algorithm was assessed with two criteria: current density and water gain rate. Verification cases produced convergence of the MEA current density { I } to within ~ 1% in a few (~ 5) iterations of the inner loop. The outer loop iterates to adjust overall hydration. Meaningful steady state, or equilibrium, is reached only when W∕ t goes below a tolerance such as 0.00005 mol m −2 s −1 . The number of iterations required varies principally from starting value of W chosen. Typically, 10-30 iterations are observed.

Model validation
The Non-equilibrium Interface Model of membrane-electrode assembly (MEA) behavior is validated through the use of a previously published experimental data set from Neyerlin et al. [19]. The interface model is intended to simulate the localized behavior of a larger device, but here is validated against a small-area PEMFC with spatially uniform flow conditions. The experiment requires well-defined MEA compositions to be meaningful. Solutions here consider only effects in the thru-plane direction and ignore variations of input parameters (properties) in the planar directions, treating those conditions as locally uniform. Figure 5 shows a flowchart of the solution scheme.
With uniform conditions within the anode and cathode gas streams, model predictions are compared against the measurements of cell ohmic resistance and the several voltage components. The voltage components are generally divided and described according to their respective loss mechanisms. Kinetic losses of the cathode's ORR reaction are presented as the first loss mechanism. Some discrepancies exist in the treatment of these losses. Ohmic losses of the membrane are the second loss mechanism. Cell ohmic resistance measurements reflect changes in hydration levels of the membrane. Ohmic losses associated with the ACL and CCL are also assessed. The measured cell voltage is shown with and without correction for high-frequency resistance (HFR).
The differential PEMFC was built with small planar dimension (i.e., 0.5 cm 2 ) and operated with gas flows of very high stoichiometry (10-100). The intention is to create aforementioned conditions as uniform as possible within the anode and cathode. The collection of experiments was performed to examine voltage losses within the cathode catalyst layer under operating conditions where full 100% catalyst utilization does not occur. The sources of voltage loss could be concurrently assessed (with catalyst layers of well-defined composition). For the purposes (here) of developing a useful interface model, a means of estimating ACL and CCL losses under all relevant operating conditions is needed. This validation step provides a means of checking the model's estimates against relevant experimental data.
The experiment explicitly made the assumption that diffusional losses within the differential PEMFC could be neglected, an assumption which sets the gas pressures, temperatures, and mole fractions from the inlet(s) of the anode and cathode as the same as those adjacent to the MEA. Diffusional losses within the catalyst layers are similarly not considered by the experiment. The task of estimating ACL/ CCL losses is then contemplated as estimating the effective protonic conduction resistance(s) within the ACL and CCL. Temperature rise was not considered. Temperature rise and diffusional losses within the MEA are estimated, however, by the interface model.

Experimental conditions
Operating conditions are shown in Table 5. Humidified hydrogen and oxygen were used as reactant gases. The gases were 100% humidified and 60% humidified at T = 80 °C. Hydrogen and oxygen partial pressures were reported as 101 kPa, and hydrogen and oxygen partial pressures were identical. The saturation pressure of water vapor is 47.79 kPa at 80 °C, and therefore, the total pressures are 149/130 kPa for the cases of 100%/60% RH. Constant gas flow rates of 1050 sccm for hydrogen and 600 sccm for oxygen feeds were used, guaranteeing minimum stoichiometric ratios of 20 (anode) and 23 (cathode) at the largest current density of 1.5 A/cm 2 . The uniformly high stoichiometric ratio of the gas flows should create spatially uniform water and current distribution within the plane of the MEA, creating the sought after differential cell flow condition. The pressure drop from inlet to outlet was reported as only 3 kPa or less; it is not considered further.
The mol fractions of hydrogen/oxygen, water vapor, and inert nitrogen are calculated from their respective partial pressures. The relevant gas pressures and mole fractions settings are shown below in Table 6.
The differential PEMFC, with active area of 5 cm 2 , was formed by two graphite interdigitated flow fields compressing the respective anode and cathode diffusion media (DM) and MEA between them. A Teflon gasket was utilized to seal the perimeter. The diffusion media were carbon fiber paper substrates (Toray, Inc.) subsequently hydrophobized with polytetrafluoroethylene (PTFE) and given an un-described surface treatment. Their thickness and composition is not further described, because that work assumes the absence of oxygen diffusion losses in these DM. Table 7 shows the MEA composition. The membrane has 1100 equivalent weight (EW) Nafion ionomer with a nominal thickness of 25 μm which is quoted at 50% RH. For modeling purposes, a dry thickness of 22 μm is used as that corresponds to a swollen membrane thickness of 25 μm at 50% RH. A crossover current density of 1 mA/cm 2 was assumed here.
The anode catalyst layer (ACL) has platinum loading of 0.35 mg Pt /cm 2 and a thickness of 12 μm. It was made from carbon-supported PT catalyst with 50% Pt/C mass ratio. The ionomer-to-carbon ratio was 1.4, which yielded nearly the given ionomer volume fraction of 0.2. The cathode catalyst layer (CCL) has platinum loading of 0.50 mg Pt /cm 2 and a thickness of 18 μm.
An estimate of the purely electronic resistances within the differential PEMFC was created experimentally. A dummy cell was created from the differential PEMFC when the MEA was removed and the device reassembled. Electrical resistance measurements of the dummy cell were used to estimate the PEMFC electrical resistances. The electronic conduction losses are considered constant, independent of temperature and humidity effects. The supplied resistance value (0.030 Ω cm 2 ) was adjusted upward, slightly, to (0.034 Ω cm 2 ) compensate for the contact resistances

Comparison of results
Interface model results are compared against the published experimental data. The experimental work reported, for the conditions given, the cell voltage V C − V A , current density I , exchange current density I x , high-frequency resistance (HFR) R , and electronic resistance R e − cnt measured from an empty "dummy cell". Diffusional resistances anywhere in the cell, whether the GDL or the cathode catalyst layer, were neglected in the work.
The HFR-corrected voltage is the sum of the measured cell voltage and measured cell ohmic losses (from HFR) as in Eq. (65). The kinetic voltage was found experimentally by adding the HFR-corrected voltage to ohmic corrections accounting for protonic conduction losses in the anode and cathode catalyst layers, in Eq. (66), where R H + ,eff ACL and R H + ,eff CCL are the effective anode and cathode catalyst layer resistances, respectively: The anode effective resistance R H + ,eff ACL was estimated from separate hydrogen-pump experiments. To determine R H + ,eff CCL , they assumed that the average membrane conductivity, measured by the high-frequency resistance measurement, could also be used to describe the average conductivity of the ionomer phase of the cathode catalyst layer. With knowledge of the thickness of each region ( t CCL and t MEM ), and the ionomer volume fraction of the CCL, CCL,Io , the effective cathode catalyst layer resistance is Eq. (67), where it was assumed in the previous work that the tortuosity of the ionomer conduction network in the CCL electrode is unity. The kinetic voltage becomes Eq. (68): The interface model was used to create results at current density values comparable to the experiment. The water gain rate in the MEA, W∕ t , is nearly zero (< 5 × 10 −5 , or several orders of magnitude lower than the water content), indicating that the water content has reached equilibrium. For the 100% RH case, W decreases with current density from anode dryness. At 60% RH, the model predicts the opposite trend. HFR measurements are compared against the interface model in Fig. 6. Experimental results at 60% RH show greater variation with current density than those at 100% RH. Though resistance readings below 0.1 A/cm 2 might be unreliable, HFR measurements at 60% RH show a drop with increasing current density. This drop is reproduced correctly by the model. It is the drop in resistance with decreasing current density, seen in the experiments, that is not reproduced by the model. It is not clear if this is due to experimental error or model deficiency. The HFR predictions at 100% RH now reflect the experimental data where resistance levels increase from 0.053 to 0.055 Ω cm 2 at 1-1.5 A/cm 2 .
Voltage measurements are compared against the interface model at 100% RH in Fig. 7. In the figure, the kinetic voltages predicted by this model utilize the active catalyst area given by the experiment. The gap between kinetic voltages and HFR-corrected voltages, predicted by this model, indicates a correct assessment of the summed ACL + CCL effective resistances. Thus, the non-equilibrium model accounts for the effective catalyst layer resistances correctly. The HFR predictions at 100% RH agreed well with the experimental values, and hence, the gap between cell voltage and the HFR-corrected voltage is also equal in both model and experiment. It follows that cell voltage levels show good agreement with measured values.
Voltage measurements are compared against the nonequilibrium interface model at 60% RH in Fig. 8. The gap between kinetic voltages and HFR-corrected voltages, predicted by this model, indicates a slight underestimation of the summed ACL + CCL losses. The HFR predictions agreed well with the experimental values, and hence, the gap between cell voltage and HFR-corrected voltage is also equal between model and experiment. The resulting cell voltage levels from the model are slightly above experimental values due to the slight underestimation of the effective ACL resistance. The differences typically were under 6-10 mV, which is similar to the amount of scatter found in the testing of the underlying cathode kinetic loss relationships [37].
Anode voltage loss, through a hydrogen-pump test, was also reported for both cases. Anode losses were reported as 30 mV at 1.5 A/cm 2 , or = 0.02 Ω cm 2 for the low-humidity case, and only reached 5 mV in the high-humidity case. Figure 9 shows anode voltage loss data against the model predictions. Anode losses of the low-humidity case are correctly assessed, while the very small anode losses of the high-humidity case were predicted as greater than the

Concluding remarks
The interface model matched experimental data produced from a well-defined differential PEMFC of known material composition. Ionic conductivity within the MEA is determined principally by the ionomer's water content. Hence, ohmic losses within the MEA are determined by the solution to a water transport problem. The model was evaluated based upon the ability to match resistance, voltage, and catalyst layer losses at varying humidity levels. With correct assessment of catalyst layer losses, it is not necessary to manipulate the treatment of cathode kinetic losses (Tafel slope and exchange current density) on a case-by-case basis to mimic the observed experimental results.
Neyerlin et al. worked to assess the kinetic behavior of the ORR (the kinetic voltage losses occurring in the cathode) [37]. They noted the converse of the problem just described. Limitations within the CCL, occurring in either the gaseous phase or as conduction losses within the ionomer, could lead the Tafel slope of the ORR to appear to increase. Earlier experiments misinterpreted the experimental results as if they indicated a change in ORR kinetics which was occurring as a (nearly) doubling of the Tafel slope of the ORR, which happened in experiments at low humidity, or where full catalyst utilization did not occur.
This interface model re-formulates the model of water transport. Electroosmotic drag and diffusion values were re-formulated according to revised water uptake curves. It also utilizes convective water transport boundary conditions which allow it to predict the hydration-dependent high-frequency resistance (HFR) as it changes with current level. Varying HFR is caused by changing MEA water contents. HFR measurements were not strictly constant with increasing current levels, indicating that membrane hydration level changes. The prior interface models incorporated Dirichlet water transport boundary conditions. They would produce a constant HFR, unchanging with current density levels.
Additional processes addressing other relevant physics were added/updated. Gas-phase diffusion and approximate thermal treatments were incorporated from recently reported research. It was also necessary to take account of ionomer thickness expansion with water uptake. Catalyst layer losses can be ascertained once the water contents, and hence ionic conductivities, are known. A reduction of ECSA, in each catalyst layer, at sub-saturated conditions, was incorporated, because it has been directly measured.
The non-equilibrium interface model thus matched experimentally derived data for high-frequency resistance (HFR), the various voltage components, and effective anode losses at reactant humidity levels of 60% RH and 100% RH. The less-important effective anode catalyst layer resistance was estimated correctly at low-humidity conditions but was slightly over-predicted for the high-humidity conditions. Good model-experiment agreement was achieved with a kinetic model whose parameters were determined from fundamental-level physical experiments reported in various other works. That is, the cathode's Tafel slope and exchange current density are not changed on a case-by-case basis; these were separately measurable. One consistent set of kinetic model parameters produced both the high-and low-humidity results.