A Physics-Based Impedance Model of Proton Exchange Membrane Fuel Cells Exhibiting Low-Frequency Inductive Loops

A physics-based impedance model of a proton exchange membrane fuel cell is developed, incorporating a coupled oxide growthoxygen reduction reaction kinetic model. The oxide layer is shown to produce a low frequency inductive loop that agrees quantitatively with the experimental inductive loop at current densities as high as 800 mA/cm2, even when kinetic and mass-transfer parameters are fit from polarization curves and cyclic voltammetry instead of electrochemical impedance spectroscopy. The importance of the inductive loop in explaining both AC and DC results is discussed. © The Author(s) 2015. Published by ECS. This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 License (CC BY, http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse of the work in any medium, provided the original work is properly cited. [DOI: 10.1149/2.0361506jes] All rights reserved.

Electrochemical impedance spectroscopy (EIS) is a widely used technique in proton exchange membrane fuel cells (PEMFCs) that separates differential contributions to cell overpotential by characteristic time constants. Analysis of EIS is difficult because there may be dozens of processes involved. It is standard to analyze EIS using equivalent electrical circuit models of the PEMFC, in which several key processes are represented by circuit elements such as the resistor, capacitor, Warburg element, and constant phase element. Fitting the experimental EIS data with an equivalent circuit model and obtaining circuit parameters is straightforward. Often, a close fit can be achieved with relatively few components, especially when constant phase elements are used. However, the challenge becomes determining which processes to include and converting the electrical circuit parameters back to meaningful cell parameters.
Although the basic building blocks of equivalent circuits, for example the Randles circuit, have an exact physical interpretation, these elements are often applied to more complex processes without accounting for the differences. The circuits may fit the data perfectly, but the parameters have lost their original meaning. Often, some transport processes are faster than double layer charging, and as a result, are inseparable from charge-transfer resistance. Additionally, slow transport processes often contribute impedance that scales with chargetransfer resistance, rather than being independent. Unfortunately, it is common in the literature to see circuit parameters reported as results without a translation to meaningful physical parameters.
An alternative approach is to model the physical processes occurring in an electrochemical cell during EIS experiments. This approach was pioneered by De Levie, 1 who calculated an analytical solution for the impedance of a porous electrode with a potential gradient, assuming linear kinetics and no concentration gradient. Keddam et al. 2 considered the alternative case of a concentration gradient without a potential gradient. Cachet and Wiart 3 considered both concentration and potential gradients for reactions following Tafel kinetics. Lasia 4 extended the De Levie model to consider Butler-Volmer kinetics, then further extended the model to consider concentration gradients in addition to potential gradients. 5 The generic porous electrode model with concentration and potential gradients has served as a foundation for device-specific impedance models.
The current state-of-the-art in PEMFC modeling was reviewed by Weber et al., 6 who emphasized the critical role of physics-based modeling in the interpretation of EIS. The subset of PEMFC models simulating impedance has been reviewed by Gomadam and Weidner 7 and more recently by Niya and Hoorfar. 8 A sampling of these PEMFC impedance models are described below. Additionally, the parallel development of physics-based impedance models in the lithium battery literature provides important context and insight, beginning with the work of Doyle et al. 9 and Meyers et al. 10 One of the first impedance models specifically for PEMFCs was developed by Springer et al. 11 They used a 1D macro-homogeneous model of the cathode catalyst layer (CCL) and gas diffusion layer (GDL) to fit a series of impedance spectra at different potentials and to demonstrate the effect of various model parameters on the modeled impedance. In the catalyst layer, gas phase and ionomer phase oxygen diffusion were lumped together into a single, 1D diffusivity. The GDL was modeled using Stefan-Maxwell multicomponent diffusion, with the assumption of saturated water vapor. The impedance was calculated by introducing a small perturbation in each variable around the steady-state solution. Assuming linearity, a new system of ODEs for the perturbation variables was produced. The impedance was calculated by solving this system of ODEs.
Springer et al. fit a set of six experimental EIS spectra over a wide range of potentials using six adjustable parameters: high frequency resistance, catalyst layer proton resistance, catalyst layer oxygen permeability, exchange-current density (geometric basis), double-layer capacitance, and GDL tortuosity. The model fit the experimental spectra well, showing two capacitive loops: catalyst-layer impedance at high frequency and GDL oxygen transport at low frequency. Single spectrum fits produced substantial variance in the model parameters, and the authors stressed that simultaneous fits over a range of potentials produce more reliable parameter estimates.
Bautista et al. 12 produced a 2D model of the MEA in order to include down-the-channel effects. The GDL was simplified by using a mass-transfer coefficient, while the catalyst layer included masstransport effects, but no ionic potential drop, following the example of Keddam et al. 2 The gas channel was modeled as a plug flow reactor. The influence of working conditions and geometric parameters on EIS spectra was investigated.
Guo and White 13 extended the standard catalyst and backing layer impedance model with a flooded agglomerate model in the catalyst layer. Previous models considered a completely flooded cathode and required unrealistically high values for oxygen permeability in the ionomer. By separating out the catalyst layer transport resistance into gas phase transport and localized ionomer phase transport, realistic values of oxygen permeability in the ionomer could be used. Additionally, by adding an extra dimension for catalyst utilization, Guo and White were able to show a quadrupled Tafel slope in cases where both agglomerate oxygen transport and ionic conduction were limiting.
Kulikovsky 14 described a simple 1D model for the catalyst layer impedance, considering oxygen diffusion, ionic conduction, and Butler-Volmer kinetics. Approximate analytical solutions were developed for limiting cases, and the numerical solution for the general case was analyzed in detail.
While the aforementioned authors have focused mainly on oxygen and proton transport, others have focused on the low-frequency behavior in EIS. It is frequently observed in the literature that the slope of the polarization curve does not match the low-frequency intercept (ca. 0.1-1 Hz) of EIS experiments. 15 Some mismatch can occur if polarization curves are done at constant gas stoichiometry, while EIS is done at constant flow rate. 16 The remaining inconsistency is explained by a low-frequency inductive loop, the beginning of which can be observed in EIS experiments extended below 0.1 Hz. 15 Several explanations for this low-frequency inductive loop have been given, including water buildup in the membrane, 17-21 buildup of ORR intermediates, [22][23][24] and platinum oxide formation from water. 23,25 The effect of water buildup in the membrane was measured experimentally by Schneider et al., 17,26,27 using a superimposed highfrequency perturbation to monitor conductivity changes during EIS. By directly measuring the changes in conductivity at the low frequencies where an inductive loop is observed, the authors conclusively proved that membrane conductivity changes were the major cause of the inductive loop for thick membranes at low relative humidity. Wiezell et al. 19 developed an MEA model further explaining this effect. All five layers of the MEA were considered, in addition to a simple flow channel mass balance. The resulting spectra included capacitive loops due to ORR kinetics and mass transport, HOR kinetics, channel oxygen depletion, and electro-osmosis-induced anode dryout. An inductive loop was observed due to the effects of water generation on conductivity and anode kinetics. The features of the simulated spectra closely matched experimental results from the same group. 20 The authors also fit their model to experimental EIS results for a range of membrane thicknesses, relative humidities, and current densities. The fitting parameters were consistent over the various conditions, and agreed with results from H 2 /H 2 operation.
While water buildup manipulates the cell Ohmic losses, ORR intermediates and the oxide layer can produce inductive loops through kinetic effects. Antoine et al., 22 Bultel et al., 24 and Roy et al. 23 modeled ORR mechanisms involving one or more intermediates and were able to simulate inductive loops. Roy et al. 23 also modeled the inductive loop through the ORR poisoning effect of platinum oxide growth. The oxide growth model was based on that of Darling and Meyers, 28 but with simpler kinetics. Mathias et al. 25 combined a twostep ORR mechanism with the Darling and Meyers model for oxide growth to further study the ORR poisoning hypothesis. The ORR mechanism involved a surface intermediate that was distinct from the slow oxide species. Characteristic frequencies for the relaxation of both the ORR intermediate and the oxide poison were derived from the model. In contrast to earlier reports, the ORR intermediate relaxation was too fast to match the low-frequency inductive loop. However, using the rate constants reported by Darling and Meyers, the oxide growth process produced an inductive loop consistent with experiments.
A better understanding of inductive loop processes is still needed. The inductive loop must be accounted for to obtain consistent parameters when fitting a model to EIS and steady-state results, especially in the case of kinetic inductive loops. As Mathias et al. 25 pointed out, in their model, the EIS and steady-state apparent transfer coefficients differ due to the relaxation of the oxide layer. If the steady-state transfer coefficient is applied to EIS analysis, kinetic resistance may be mistaken for mass-transport resistance. Additionally, model validation is under-addressed in the literature, with only a handful of studies making quantitative comparisons between impedance models and experimental results. Most quantitative comparisons are merely model fits, and not separate validation experiments. Fitting without validation increases the likelihood that EIS losses are assigned to the wrong process due to a missing effect in a model.
The present work builds on earlier models, combining the usual porous electrode transport phenomena with oxide poisoning, heat generation, gas channel pressure drop, and gas stoichiometery. Model parameters are sourced from independent experiments, and the model is validated over a range of conditions. Water buildup and the oxide layer are both found to be capable of producing inductive loops, but only the oxide layer makes a significant contribution under the conditions tested. Additionally, the cell hardware is shown to influence EIS spectra through flow channel pressure drop, flow stoichiometry, and flow field heat transfer. This paper describes the PEMFC model, the procedure for fitting parameters, and EIS results. An overview of the inductive loop processes is also provided. A more detailed discussion of the inductive loop processes will be provided in a second publication, with an exploration of the parameters and conditions that determine the magnitude and characteristic frequencies of the inductive impedance.

Model
The model consists of eight parts representing the five layers of the membrane electrode assembly (MEA), the two flow fields, and a frequency response analyzer (FRA). The model was implemented in gPROMS v3.5.3 (Process Systems Enterprise Ltd.). For brevity, a description of the model physics without equations is given in this section, and the detailed equations are given in the Appendix.
The MEA model is a 1D continuum model, with Stefan-Maxwell diffusion in the gas phase, concentrated solution theory in the ionomer phase, and Ohm's law in the carbon phase. Ionomer transport parameters are a function of ionomer water content. The anode catalyst layer is modeled as a fully reversible, planar electrode. The CCL includes an additional dimension of oxygen transport through the flooded agglomerate with thin film model, which utilizes the pseudo-steady-state approximation. The oxygen reduction reaction (ORR) follows Tafel kinetics with an additional oxide layer effect, as described later in this section. As the model is 1D, the ORR rate changes with overpotential, oxygen concentration, and oxide coverage through the thickness of the catalyst layer. Oxide growth and double-layer capacitance also contribute to the current. Convection in the GDL is not modeled explicitly, but its effect is approximated by reducing the tortuosity parameter to match the mass-transfer resistance. Knudsen diffusion is neglected. While Knudsen diffusion may account for a significant fraction of the gas-phase mass-transport resistance in the catalyst layer, 29 the overall effect would be small because the total gas-phase mass-transport resistance is small in the catalyst layer. Water condensation is not included in the model, and care is taken to restrict the present analysis to sub-saturated conditions. Joule heating, membrane hydration, ORR overpotential, and halfreaction enthalpies are all accounted for as heat-generation sources, and heat transfer occurs through conduction. The heat capacity of each layer is included, although the transient effect is negligible except for the flow fields. The most significant effect of temperature is on relative humidity, resulting in ionomer dryout. All other effects of temperature are assumed negligible, including effects on transport and kinetic parameters.
The flow fields are modeled with separate 1D domains for channel flow and heat transport. Channel flow is modeled with a 1D mass balance, assuming a uniform gas flux into the GDL determined by the 1D MEA model, which lacks a down-the-channel dimension. Flow channel pressure drop is accounted for by an empirical correlation, and a linear pressure profile is assumed. Partial pressure boundary conditions for the MEA model are determined from the average over the 1D channel at any given point in time. This scheme approximates the important 2D effects without the additional model complexity. The 1D heat transport domain is needed to account for the thermal resistance of the thick graphite flow fields in typical research hardware. An additional heat-transfer resistance is applied at the outer boundary to account for the electrical insulator between the current collectors and the temperature-controlled end plates. Unlike in the thin MEA, the heat-transfer time constant in the thick graphite flow field is slow enough to be relevant in EIS.
The platinum oxide layer is critical to the model results. Following the work of Redmond et al., 30 the oxide model includes a fast, chemisorbed oxide and two slow, place-exchanged oxides (planar-site and edge-site oxide). The slow oxides block chemisorption sites, such that the initial reversible chemisorbed oxide is gradually replaced with irreversible place-exchanged oxide. In order to account for the oxide F521 layer effect on ORR kinetics, the ORR rate equation is first order in vacant chemisorption sites. Other researchers have used a double-trap model for ORR kinetics and oxide growth, 31,32 in which oxygen is reduced through two parallel routes involving adsorbed intermediates. The oxide layer consists of these adsorbed intermediates, and their buildup at high potentials causes the dual Tafel slope. However, in experiments, oxide growth continues over very long timescales, 33 which are not compatible with the short timescales for ORR intermediates required to support a reasonable reaction rate. 25 For this reason, oxide growth was treated as a separate process and modeled independently of ORR using the Redmond et al. model.
Impedance is calculated following the approach of Boaventura et al. 34 EIS is simulated in the time domain, and the results are transformed to the frequency domain following the operating principles of an analog frequency response analyzer (FRA). 21 The FRA uses a finite Fourier transform, to determine the magnitude and phase of the response at a given frequency according to the following process. Real and imaginary excitation signals which are 90 • out of phase with each other are defined. The cell voltage is multiplied by these excitation signals and integrated using an ODE built into the model. The same procedure is done for the cell current, and the ratio of the complex voltage and current determines the impedance. The time-domain approach allows easy implementation of EIS with any transient model.
The model parameters which are sourced from the literature or from manufacturer data sheets are listed in Table I. The remaining parameters are fit by experiment as described in the Results and Discussion section.

Experimental
For parameter fitting and model validation, experiments were performed on a commercial 25 cm 2 MEA from Ion Power, Inc. The MEA used a Nafion 212 membrane and a Pt/C catalyst with 0.3 mg Pt /cm 2 loading on each electrode. Toray carbon paper (TGP-H-060, 0.19 mm thickness) with 5% PTFE wet-proofing (Fuel Cell Store) was used as a GDL. A 0.16 mm PTFE gasket was used to set the GDL compression. The cell was tested in single cell hardware from Fuel Cell Technologies with Poco graphite flow fields in a triple serpentine flow pattern. The approximate dimensions were 0.7 mm channel width, 0.9 mm land width, and 1.1 mm channel depth. The flow fields were backed by gold-plated current collectors, which were separated from the temperature-controlled aluminum end plates by a PTFE-impregnated fiberglass tape (SG25-05, Saint-Gobain Performance Plastics, 0.19 mm thickness). The cell was controlled using a Scribner Associates 850e fuel-cell test system with a Model 880 frequency response analyzer. Unless otherwise noted, all experiments were conducted at 80 • C with inlet relative humidity of 75%. No backpressure was applied.
Three polarization curves were performed using different cathode gases: 1% oxygen in nitrogen, 1% oxygen in helium, and air. Polarization curves consisted of a series of 15 minute constant current holds, starting with the highest current first. The potential was averaged over the last five minutes of each hold. After each polarization curve, galvanostatic EIS was performed at a range of current densities. Each EIS sweep was preceded by a brief hold at a low potential and a 15 minute hold at the DC current density. The AC amplitude (RMS) was 5% of the DC current density, and the frequency was swept logarithmically from 10 kHz to 10 mHz with 10 points per decade. Under air, the polarization curve and EIS used anode and cathode flow rates of 1.49 mmol/s (2000 sccm) and 3.72 mmol/s (5000 sccm), respectively. The experiments conducted with 1% oxygen used anode and cathode flow rates of 0.372 mmol/s and 3.72 mmol/s, respectively.
A Tafel plot was obtained by cyclic voltammetry (CV) at 0.5 mV/s from 0.7 V to 0.95 V with oxygen as the cathode gas. The anode and cathode flow rates were 1.49 mmol/s and 3.72 mmol/s, respectively. The high frequency resistance (3 kHz) was used to determine the iRcorrected potential. Potentiostatic EIS with nitrogen at the cathode was used to measure membrane and ionomer conductivity. The potential was 0.4 V with an amplitude of 5 mV RMS , and the hydrogen and nitrogen flow rates were 74 μmol/s. The frequency range was 10 kHz to 100 Hz with 20 points per decade, and 100 Hz to 1 Hz with 10 points per decade. The inlet relative humidity was varied over a range from 42% to 100%, and a total of 28 EIS spectra were recorded. Using Matlab R2010b, the catalyst layer proton resistance was fit to a transmission line model, with five fitted parameters: cable inductance, L wire , membrane resistance, R mem , catalyst layer resistance, R cl , and double layer constant phase element (CPE) admittance and exponent, Q dl and φ. The derivation of Equation 2 is provided in the electronic supplementary information (ESI) in the section "Derivation of H 2 /N 2 impedance." Two CVs were conducted to measure electrochemically active surface area (ECA) and oxide growth parameters. The CVs were performed with hydrogen and nitrogen flow rates of 74 μmol/s and 37 μmol/s, respectively. ECA was determined from a CV at 25 • C and 100% RH, using a sweep from 0.05-0.60 V at 50 mV/s. The cathodic hydrogen adsorption peak was integrated between the baseline near 0.4 V and the local maximum near 0.09 V, assuming a charge of 210 μC/cm 2 . 35 The oxide growth parameters were determined from a CV at 80 • C and 75% RH, using a sweep from 0.05-1.00 V at 50 mV/s.

Results and Discussion
A number of approaches to parameter fitting and model validation are possible. Here, several experiments are used that are designed to isolate the most important model parameters, and the model is validated with EIS spectra over a range of current densities and oxygen concentrations. The fitting experiments are: H 2 /N 2 EIS for conductivity, CV for ECA and oxide parameters, Tafel curve for ORR rate constant (from slow CV under oxygen), and limiting current density analysis for mass-transfer parameters (using polarization curves under 1% oxygen in nitrogen and helium). It is important to stress that the validation experiments, EIS with varying oxygen concentrations, were not used to fit any model parameters.
Conductivity.-To begin, ionomer conductivity and catalyst-layer tortuosity are fit using EIS measurements with nitrogen at the cathode. 15,36,37 A total of 28 EIS spectra were acquired at nine values of relative humidity ranging from 42% to 100%. Figure 1a shows representative experimental EIS spectra with the modeled spectra overlaid. The Nyquist plot shows the typical shape for a transmission line with no Faradaic process: a 45 • line at high frequencies transitioning into a (nearly) vertical line at low frequencies. The high frequency intercept is the membrane resistance plus electronic resistances and contact resistances, which are assumed to be negligible. If the low frequency line is extrapolated down to the x-axis, this intercept has been shown to be equal to the high frequency resistance plus one-third of the catalyst layer ionomer resistance. 15 Although the graphical reading is easiest to understand and is a correct interpretation of the transmission line model, a more precise procedure is to fit the entire spectrum with the transmission line model of Equation 2.
Each of the 28 experimental spectra was fit using Equation 2, and the resulting membrane conductivity, effective catalyst layer conductivity, and catalyst layer ionomer tortuosity are shown in Figure 1b. The swelling effect was ignored, and the thicknesses used in the calculation were 50 μm for the membrane and 17 μm for the catalyst layer. The remaining fit parameters are not used in the PEMFC model, but are listed in Table S1 of the ESI for reference. Except at low relative humidity, the CPE phase parameter was ca. 0.98, indicating nearly ideal double-layer capacitive behavior at 0.4 V on Pt/C. Catalyst layer tortuosity was determined from the ratio of the membrane and catalyst layer conductivities, after accounting for the ionomer volume fraction of 0.16. The ionomer volume fraction was calculated from the measured porosity of 0.65 assuming an ionomer to carbon ratio of 1:1, as detailed in the ESI. The catalyst layer effective conductivity exhibits stronger water activity dependence than the membrane conductivity, which is consistent with literature reports on conductivity, 36 but at first glance appears contradictory with literature reports finding lower ionomer hydration in catalyst layers. 38 However, even with reduced hydration levels, the catalyst layer ionomer will swell with increasing relative humidity. The resulting increase in volume fraction and decrease in tortuosity appears to dominate any difference in bulk vs. thin-film ionomer conductivity. At high water activity, the calculated tortuosity falls below one, although the true tortuosity after accounting for swelling would be larger than one. The data were fit with an exponential function to obtain the correlations used in the PEMFC model, which are listed in Table II. The data at 100% RH were excluded to achieve a better fit at subsaturated conditions.
Kinetic parameters.-Cyclic voltammetry was used to determine ECA and the growth of the oxide layer on platinum. ECA was calculated using the charge for hydrogen adsorption, assuming a charge of 210 μC/cm 235 The room temperature CV is included in the ESI, Figure S1, from which an ECA of 79 m 2 /g was obtained. Oxide growth   Figure 1).

Symbol
Name Value is detailed in Figure 2, which consists of experimental and simulated cyclic voltammograms at 80 • C, 75% RH, and 50 mV/s scan rate.
The experimental curve has been shifted by the crossover current, 2.5 mA/cm 2 . The oxide growth parameters were fit from the CV and are listed in Table III. Figure 3 shows experimental and simulated Tafel plots with oxygen at the cathode. The Tafel plot was acquired through a CV from 0.70  Frumkin interaction energy edge oxide 550 kJ/mol to 0.95 V at 0.5 mV/s, with a low potential prehold to reduce the oxide layer. The experimental Tafel plot is corrected for hydrogen crossover of 2.5 mA/cm 2 and high frequency resistance, measured at 3 kHz. The simulated curve is corrected for membrane resistance.   Figure 2 reproduce the hysteresis in the Tafel plot nearly perfectly. Additionally, the simulated Tafel slope on an oxide free surface is 140 mV/decade, but the Tafel slope is reduced to ca. 70 mV/decade by the oxide interactions, matching the experimental Tafel plot. The transition in Tafel slope matches the experimental observation of a dual Tafel slope 39 on polycrystalline platinum, although the transition occurs at too low of a potential to be seen in Figure 3, which is consistent with reports on MEAs. 40 Mass-transfer parameters.-Limiting current density analysis was used to calculate the mass-transfer parameters. The three masstransfer parameters that require fitting are agglomerate radius, thin film thickness, and GDL tortuosity. Although literature values are available for GDL tortuosity, 41 the model uses an adjusted parameter which accounts for GDL convection as well. The adjusted GDL tortuosity is the value which produces a mass-transfer resistance through diffusion alone that is equivalent to the actual mass-transfer resistance due to the combined effects of diffusion and convection. This method assumes an equal distribution of mass-transfer resistance through the GDL, which is not strictly valid because convection and diffusion have separate driving forces. Improved convection models would lead to different concentration profiles in the GDL. In reality, the convection effect should be stronger near the channel than near the catalyst layer, essentially changing more than D ef f . As a result the diffusion time constant ( 2 /D ef f ) may be overestimated. However, in these simulations, the effect on impedance is minimal because GDL mass transfer is already too fast to resolve from the charge-transfer semicircle. Flooded agglomerate parameters are specific to the catalyst layer fabrication technique and also require fitting. Limiting current density analysis is capable of separating ionomer and gas phase mass-transfer resistance. 42,43 However, to separate the two ionomer mass-transfer parameters, agglomerate radius and film thickness, an additional datum is required: the cell potential at 100 mA/cm 2 .
For low oxygen concentrations, the mass-transfer limiting current density can be reached with minimal ohmic losses and without significant water buildup or the risk of condensation. Due to the cell hardware, differential conditions were not obtainable; the oxygen stoichiometry was as low as 2.6. Therefore, the down-the-channel concentration gradients need to be accounted for. A simplified 1D model (down the channel) is used in which the gas phase and ionomer phase are treated as mass-transfer resistances, and the oxygen partial pressure at the catalyst surface is zero. The details of the model are provided in the ESI, Figure S3 and associated text. Accounting for the pressure drop in the gas flow channels, the outlet mole fraction is ln where k p,i and k x,g are the ionomer and GDL mass-transfer coefficients, N ch is the molar flow rate per unit MEA area, and P(y) is the flow channel pressure at point y, the normalized distance from the inlet. While k p,i is expressed in terms of a partial pressure driving force, k x,g is expressed in terms of a mole fraction driving force, as gas phase diffusivity is inversely proportional to pressure. Assuming a linear pressure profile, Equation 3 becomes [4] The gas phase and ionomer phase mass-transfer coefficients may be separated by comparing the limiting current density for O 2 /N 2 and O 2 /He mixtures. For the O 2 /He mixture, k x,g is replaced by α D k x,g , where α D is the ratio of oxygen diffusivity in humidified helium to diffusivity in humidified nitrogen. The oxygen diffusivity in a gas mixture is estimated from the binary diffusivities according to The limiting current density was measured at 75% RH with a dry gas flow rate of 3.72 mmol/s at the cathode and 0.372 mmol/s at the anode. The observed pressure drop was 71.6 kPa for nitrogen and 47 kPa for helium. Polarization curves under these conditions are shown in Figure 4. Somewhat surprisingly, the limiting current density was 223 mA/cm 2 for nitrogen, and only 214 mA/cm 2 for helium, indicating that the higher pressure in the nitrogen case outweighs the lower gas-phase diffusivity. The calculated mass-transfer coefficients were k p,i = 8.77 × 10 −3 mol/(m 2 · s · kPa) and k x,g = 7.55 mol/(m 2 · s) (humidified nitrogen). Ignoring convection, the tortuosity of the GDL is determined from k x,g by where D O 2 is the diffusion coefficient of oxygen in the humidified nitrogen mixture and G is the thickness of the GDL. The calculated tortuosity was 0.62, compared to a typical literature value of 2.85 41 While τ = 0.62 is suitable for the present model to account for convection, it is not a true measure of the tortuosity of the GDL. For an agglomerate of radius, R agg , including an ionomer film of thickness d f , the ionomer phase mass-transfer coefficient is The agglomerate radius was fit by matching the potential at 100 mA/cm 2 for the O 2 /N 2 case to the model, with the film thickness determined from Equation 7. The fitted agglomerate radius and film thickness are 420 nm and 27.1 nm, respectively. The experimentally determined model parameters are listed in Table II. The modeled agglomerate size is large and would be easily observable by SEM. However, such large agglomerates are never observed. 43 Such overestimates of agglomerate dimensions are common in the literature, 13,[44][45][46] and may result from assuming bulk permeability values for ionomer thin films. By extrapolating measurements of masstransport resistance to zero thickness in thin ionomer films, Suzuki et al. 43 concluded that an interfacial mass-transfer resistance is present. When including this extra transport resistance, the modeled agglomerate dimensions were roughly the size of the primary carbon particles,  and consistent with catalyst layer images. Furthermore, the present model required the unrealistic assumption that agglomerate tortuosity was unity in order to achieve the correct balance between external and internal mass-transfer limitations. In the present model, it is better to treat the agglomerate model as a purely empirical description of internal and external mass-transfer limitations, and not to treat it as a literal model of the catalyst layer microstructure.
Model validation.-With the model parameters determined from independent experiments, the model is evaluated by comparing simulated and experimental EIS spectra under a variety of conditions. Figure 5 shows EIS results for a cell running on 1% O 2 /N 2 . The cell temperature was 80 • C, and the dry gas flow rates were 0.372 mmol/s anode and 3.72 mmol/s cathode, both at 75% relative humidity. In addition to simulated and experimental EIS curves, the DC resistance from the polarization curve is included for comparison as filled symbols on the real axis labeled "DC". Figure 5a displays the Nyquist plot for low current densities, and Figure 5b shows the Nyquist plot for high current densities. Figure 5c and Figure 5d are low and high current density Bode plots, respectively, showing the imaginary component of impedance. The spectra were split between "low" and "high" current densities to avoid overcrowding of the plots, but no physical significance in terms of a characteristic current density for a particular process was intended. Each spectrum consists mainly of a single capacitive loop, stretching from the high frequencies to the traditional low frequency intercept near 0.1 Hz. However, there is a significant difference between the ca. 0.1 Hz intercept and the DC resistance. At the higher current densities, the clear beginning of a low frequency inductive loop is observed. Thus, the 0.1 Hz intercept will be referred to as the intermediate frequency intercept to distinguish it from the true low frequency intercept, which is hard to observe in EIS, but can be measured from the polarization curve. Presumably, if it were practical to extend EIS into the μHz range, one would observe the same low frequency intercept at the end of the inductive loop as is measured from the polarization curve.
The modeled spectra show matching features to the experimental results, with a capacitive loop followed by an inductive loop. At high frequencies, the model matches the length and position of the 45 • line caused by porous electrode effects and even exhibits non-ideal capacitive behavior. In the model, the non-ideal capacitive behavior, usually requiring an equivalent circuit fit with a CPE, is caused by growth of the oxide layer, for which heterogeneity produces a dispersion of time constants. For current densities up to 50 mA/cm 2 (Figure 5a (Figure 5c,5d) show that the model underestimates the characteristic frequencies by a factor of 2-3. The frequency mismatch is due to extra model pseudocapacitance from the chemisorbed oxide, which must be somewhat overestimated by the oxide growth model for long hold times.
Under hydrogen and air, the EIS spectra are shown in Figure 6. The cell temperature was 80 • C, and the dry gas flow rates were 1.49 mmol/s anode and 3.72 mmol/s cathode, both at 75% relative humidity. Similar to the 1% O 2 case, the current densities up to 400 mA/cm 2 show excellent agreement in the Nyquist plot (Figure 6a), but underestimate the frequencies in the Bode plot (Figure 6c). The results at high current density (Figure 6b) are mixed, with good agreement at 800 mA/cm 2 . At 1600 mA/cm 2 , the model underestimates the impedance. The Bode plot (Figure 6d) shows that the frequency mismatch begins to be resolved at higher current densities. The model does underestimate the high frequency resistance by about 15%. Whether this is due to RH droop in the humidifiers at maximum flow rate or a model deficiency is unknown. Interestingly, at 1600 mA/cm 2 , the model shows a small peak at low frequencies on the Bode plot. This peak is caused by heat generation and storage in the flow field blocks.
The inductive loops observed in Figure 5 and Figure 6 are caused by the relaxation of the oxide layer on platinum. An additional source of inductive loops in PEMFCs is water generation and storage in the membrane. 17,[19][20][21]26 An isothermal model predicts an inductive loop due to product water buildup. As current density increases, the increased water production causes a drop in the membrane and ionomer resistance, but only if the frequency is low enough for the membrane to absorb the additional water. This modulation of the membrane resistance is effectively negative impedance, thus causing an inductive loop. However, when heat transfer is considered, the cell temperature rise at higher current densities partially or fully offsets the increased water generation. The hydration response of the membrane may be reversed and result in a capacitive loop and positive impedance. In the present model, heat generation offsets the water buildup to the extent that a small capacitive feature is observed in the model at 1600 mA/cm 2 . In contrast, the experimental measurements of Schneider et al. 17,26 were performed at lower temperature, lower inlet relative humidity, lower flow rates, and with a larger cell, resulting in a large inductive loop from water generation. Thus, while water generation is not responsible for an inductive loop under the conditions in the present work, it can be very important under other conditions. The results of Figure 5 and Figure 6 are summarized in Figure 7 and Figure 8, which show the intermediate frequency and DC inter-cepts at each current density. As mentioned previously, both intercepts are matched almost perfectly up to 25 mA/cm 2 in 1% O 2 /N 2 , and up to 800 mA/cm 2 in air. At higher current densities, the inductive loop is too small to fit both intercepts. At high current densities, one might suspect that localized condensation could occur and add to the transport losses. However, even at 2 A/cm 2 , the model predicts only a   2.6 kPa gradient in water partial pressure from the catalyst layer to the flow channel. The closest the channel partial pressure comes to the vapor pressure is 13 kPa. While there will be spatial variations in the GDL mass-transfer resistance due to the serpentine flow pattern, the margin should be sufficient to prevent any localized condensation under these conditions. Finally, the polarization curve under air is shown in Figure 9. The steady-state polarization curve is not a challenging test considering that the mass-transfer parameters were fit from a steady-state polarization curve, albeit at 21 times lower oxygen concentration. However, the air polarization curve does demonstrate that the model performs well under both DC and AC conditions.

Conclusions
In this work, a physics-based PEMFC impedance model was demonstrated with additional physics to account for the oxide layer, heat generation, and cell hardware effects. With these extra effects, the model matches more of the trends and features of experimental EIS measurements, even while using very different experiments for parameter fitting and model validation. The proposed ORR model, which links the previously published oxide growth model 30 to ORR kinetics through vacant chemisorption sites, is shown to produce a very low frequency inductive loop. This inductive loop quantitatively agrees with the observed inductive loop at current densities up to 800 mA/cm 2 under air and 25 mA/cm 2 under 1% O 2 /N 2 . The model also can generate an inductive loop through the effect of water buildup, although in the cases studied, heat generation completely offsets the water buildup.
The largest impediment to accurate physics-based impedance models has been the low frequency inductive loop. Without accounting for the processes responsible for this feature, a model can only fit EIS spectra (>0.1 Hz) or steady-state experiments, but not both. The present model was able to explain most of the observed inductive loop, as well as the ORR Tafel slope, ORR hysteresis, and CPE behavior. These results indicate that the oxide layer has a major effect on PEMFC impedance and must not be ignored in modeling.

Acknowledgments
Toyota Motor Engineering & Manufacturing North America supported this work.

Appendix: Detailed Model Specification
The CCL is the most complex sub-model, and consists of mass, charge, and heat balances, transport equations, and kinetic equations. The GDLs and membrane are modeled with a subset of the CCL model, excluding the absent phases. The CCL model is described in detail below.
The conserved variables are gas partial pressures, concentration of water in the ionomer, double layer charge, and temperature. The seven balance equations for oxygen, water vapor, nitrogen, ionomer water, electronic charge, ionic charge, and heat are: Several key symbols will be defined as they are introduced, but all symbols are included in the List of Symbols section. The concentration of water in the membrane is expressed as moles of water per equivalent of acid, λ. The generation term, j ORR , is the rate of the oxygen reduction reaction (ORR) on a volumetric basis, while j v , is the rate of water evaporation from the ionomer. Different transport equations are used for each phase. In the gas phase, the Stefan-Maxwell equations, describe multicomponent diffusion. D ij is the binary diffusion coefficient for a pair of gases and is estimated from kinetic theory 47 as Pressure, P, is measured in bar, and average molecular diameter, σ ij , is measured in nm to give D ij in units of m 2 /s. Parameters were obtained from Ref. 48. The Stefan-Maxwell equations, [A8], and the requirement that total pressure is the sum of all partial pressures together yield three independent equations.
The ionomer phase is described using concentrated solution theory, as detailed in Fuller and Newman. 49 The three transport properties are conductivity, κ, water diffusivity, α, and electro-osmotic drag coefficient, ξ. The transport equations are The subscript, eff, denotes effective transport properties, which take into account the porosity and tortuosity of the relevant phase. The chemical potential referenced diffusion coefficient, α, was determined from the Fick's law diffusion coefficient of water, D w,i , according to [A12] The gas species are described by the ideal gas equation of state, and water in the ionomer is described by the empirical correlation of Springer et al., 50 The double layer charge is determined by the interfacial potential difference, The ORR rate equation, is based on the Tafel equation, with a first order dependence on vacant catalytic sites, v , and a power-law dependence on oxygen partial pressure. In order to model the effect of the oxide layer, it is necessary to model the growth of the oxide layer. A previously developed oxide-growth model 30 is used with some minor adjustments. Although this is not the first time vacant sites have been used to explain the oxide layer ORR hindrance, 39 most previous attempts have used simple equilibrium oxide coverage models, which do not match experimental oxide coverage results very well. The oxide growth model of Redmond et al. matches growth rates and coverages over a wide range of timescales.
The oxide layer is formed according to the reaction scheme which consists of three components: chemisorbed OH, place-exchanged PtO 2 on planar sites, and place-exchanged PtO 2 on edge sites. A brief introduction to the notation will be helpful. Surface concentration (mol/m 2 ) is represented by , where 0 is a constant equal to the surface concentration of planar-site platinum atoms (excludes edges), and where X e is the fraction of platinum surface atoms on edges. Both PtO 2 species exhibit heterogeneity in formation energy represented by μ. The coverage at each formation energy is x 2 or x 3 , related to total coverage by, The degree of heterogeneity is determined by the initial distribution of oxides formed in the anodic process, ψ. A normal distribution with variance σ 2 and σ 3 is used for planar and edge oxides: [A25] Both forms of PtO 2 disrupt chemisorption sites according to (cf. eq. 10 in Ref. 30 where χ is the number of sites disrupted by a single unit of PtO 2 . In Ref. 30, only planar-site PtO 2 is assumed to disrupt chemisorption sites. Equation A26 applies given the assumption that PtO 2 units are placed randomly on the surface and may overlap. Chemisorption is a fast, quasi-equilibrated reaction and follows the equilibrium expression (cf. eq. 15 in Ref. 30), The rate equation (cf. eq. 16 in Ref. 30) for planar PtO 2 is A small change in the anodic Frumkin term has been made from Ref. 30, with θ 2 , the total coverage, being replaced by , the coverage at a particular μ relative to the distribution function. During the course of normal CV simulations, this change is inconsequential; however, after long holds the former version results in a very inert oxide layer, while the latter version retains some reactivity.
Similarly, for the edge site oxide, the rate equation is (cf. eq. 18 in Ref. 30) The current due to oxide formation is Finally, the surface concentration of vacant chemisorption sites, required in the ORR rate equation A15, is v = (1 − θ 1 ) * . [A31] To model oxygen transport through the ionomer in the catalyst layer, a flooded agglomerate model with an extra ionomer film is used. Assuming an m th -order reaction, an approximate analytical solution is possible. First, an equivalent homogenous rate constant, k is defined as Next, the Thiele modulus, φ, is defined for the active core of the agglomerate as where p O 2 , f is the partial pressure of oxygen at the film-agglomerate core interface, and D O 2 ,i,ef f is the effective diffusivity of the agglomerate core. The effectiveness factor, which is the ratio of the average rate to the rate in the absence of mass-transfer limitations, is a function of the Thiele modulus, The film transport resistance causes a drop in p O 2 , f , according to With the effectiveness factor, rate constant, and film partial pressure defined, the macroscopic homogeneous reaction rate is calculated as [A36] The transport properties of the ionomer are dependent on the level of hydration. Conductivity was determined experimentally as detailed in the Results and Discussion section. The diffusion coefficient of water is taken from Fuller and Newman, 49 The electro-osmotic drag coefficient was determined previously by the streaming potential method, 51 . [A38] ) unless CC License in place (see abstract In addition, several test hardware effects are considered. The changing composition of the gas in the flow channel is considered using a 1D mass balance, where y is the dimensionless distance down the flow channel, ch ef f is the ratio of flow channel volume to MEA area, N ch,i is the molar flow rate of gas i in the channel, normalized by MEA area, and N G i is the flux of gas i into the flow channel from the GDL. As the MEA model is 1D, a uniform gas flux into the GDL must be assumed down the length of the channel. A linear pressure drop profile is assumed down the length of the flow channel, and the inlet total pressure is determined by an experimentally measured correlation, [A40] using the parameters listed in Table AI. The measured pressure drop is shown in the ESI, Figure S2. The partial pressures and molar flow rates are related by . [A41] The GDL boundary partial pressures are set to the average partial pressures in the channel. This boundary condition couples the time-varying channel composition to the MEA model, capturing the dynamic effects of finite gas stoichiometry. Two heat-transfer resistances in the test hardware are considered. The graphite flow field is modeled with a 1D heat balance, where A c P and A k are the effective cross-sectional areas for heat storage and heat conduction, accounting for the larger cross-sectional area of the flow field than the MEA (58.1 cm 2 vs. 25 cm 2 ). The thin, electrically insulating sheet between the current collector and the end plate also adds a small thermal resistance, which is modeled as a boundary condition rather than a 1D domain. The boundary condition is where h ins,eff is the effective heat transfer coefficient of the insulating layer, and T 0 is the regulated temperature of the end plates. The effective cross-sectional areas were determined with a steady-state 3D model of heat transfer in the flow field and electrical insulator. The ratio of heat fluxes and heat storage in the 3D model to the 1D model were used to calculate the effective cross sectional areas listed in Table I.