One-Dimensional Impedance of the Cathode Side of a PEM Fuel Cell: Exact Analytical Solution

A transient one-dimensional model for the coupled performance of the cathode catalyst layer (CCL) and the gas-diffusion layer (GDL) in a PEM fuel cell is linearized and an exact analytical solution for the cathode side impedance Z is derived. The solution is valid at a low current density, when the proton and oxygen transport losses in the CCL are not large. The zero-frequency resistance contains the faradaic term, the terms due to the proton and oxygen transport in the CCL and due to the oxygen transport in the GDL. The analytical equation for Z is fast enough for least-squares ﬁtting of experimental spectra, which yields the GDL oxygen diffusivity, the Tafel slope of the oxygen reduction reaction, and the CCL proton conductivity and oxygen diffusivity. As an example, the model is ﬁtted to three measured spectra of a high-temperature PEMFC; the ﬁtting parameters are discussed. an the terms of the Commons

Polymer electrolyte membrane fuel cells (PEMFCs) are close to commercialization: leading car manufacturers claim beginning of mass production of fuel cell cars within the next few years. Yet, a lot of research has to be done to improve PEMFC durability and performance.
One of the most powerful in situ techniques for fuel cell characterisation is electrochemical impedance spectroscopy (EIS). 1 In a working fuel cell, the transport and kinetic processes affect the cell impedance, as all these processes are linked to the charges production and consumption in electrochemical reactions. A great advantage of EIS is that it nicely separates potential losses corresponding to the processes with different time scales, whereas the cell polarization curve shows only the sum of these losses. Interpretation of EIS spectra requires modeling.
Historically, EIS has been developed for electrochemical studies on a plane single-crystal electrode surface immersed in a liquid electrolyte. 1,2 A number of classic solutions for the impedance of this system in different conditions has been derived e.g., Warburg transport impedance, Gerischer impedance for the system with transport and chemical reaction in the electrolyte 3 etc. The Gerischer impedance has been used for analysis of impedance spectra of SOFC anode, in which the electrochemical conversion occurs at the anode/electrolyte interface. 4,5 However, PEMFC electrode is a porous system with the reaction distributed over the electrode depth. It can be shown that the impedance of this system is not equivalent to the Gerischer impedance of a planar electrode (see below).
Analytical solutions for impedance of a single cylindrical pore have been derived by Lasia. 6,7 DeLevie suggested a transmission line model for the porous electrode and discussed a physical model for the electrode impedance. 8 However, the transmission line of deLevie does not include oxygen transport elements and equivalence of this simple model to the physical model resulting from the conservation equations has not been rigourously proofed. Every physical impedance model is a tool for fitting experimental impedance spectra in order to obtain physical parameters of the system. Accurate least-squares fitting algorithms involve numerous calls to subroutines for calculation of the real and imaginary part of impedance. These subroutines must be fast, otherwise the code would be time-consuming and not suitable for massive spectra fitting. This limits the use of numerical impedance models in fitting algorithms: as soon as the code includes numerical solution to a boundary-value problem for the overpotential and/or the oxygen concentration in the CCL, it is slow. Particularly time-consuming is calculation of impedance at high frequencies, when the perturbation is "trapped" close to the membrane interface and it exhibits very high spatial gradients. As a result, in spite of development of numerical physical models for the CCL impedance, [9][10][11][12][13][14][15] fuel cell developers routinely use the standard library of transmission line elements inherited from the classic EIS studies for impedance spectra fitting.
In recent years, following a milestone paper by Eikerling and Kornyshev, 16 a number of "fast" analytical solutions for impedance of the porous catalyst layer with polymer electrolyte have been derived, [17][18][19][20][21][22][23] (many useful references on the PEMFC impedance studies can be found in Ref. 24). So far, however, no complete analytical impedance model for the PEMFC cathode side, which would include impedance of the GDL has been developed.
The model of Springer et al. 25 combines a numerical model for the CCL impedance and an analytical solution for the GDL impedance. However, the equation for the GDL impedance is derived in Ref. 25 neglecting the oxygen transport in the CCL, as if the GDL were connected to the infinitely thin catalyst layer with no transport losses. The limits of validity of this approximation are not clear. Accurate approach requires to solve the CCL and GDL impedance problems simultaneously, taking into account continuity of the oxygen concentration and flux at the CCL/GDL interface (see below). In addition, due to the numerical sub-model for the CCL, the impedance model of Ref. 25 is hardly suitable for massive spectra fitting.
A combination of analytical and numerical approaches to calculation of the GDL impedance has been developed by Bultel et al. 11 The model is based on the Stefan-Maxwell transport equations for the air components. The oxygen mass balance equation in the GDL takes into account the effect of Tafel slope doubling at high cell currents due to the proton transport losses in the CCL. However, the CCL impedance itself is calculated in a low-current approximation, under the assumption of constant through-plane shape of the local overpotential. Nonetheless, at real cell operating currents, the model 11 generates impedance spectra which are similar to the experimental spectra. It is worth mentioning that no attempts to fit measured spectra have been done in this work.
Below, we develop a physical model for the impedance of the system "CCL + GDL" in a PEM fuel cell. The model is an extension of a recent model for the CCL impedance. 22 In the case of a small cell current, the model is solved analytically, which allows us to construct a fast code for impedance spectra fitting. Examples of such fitting are given and the set of cell parameters resulting from fitting is discussed. The Maple worksheet with this code is available for download; the link is indicated in the Results and Discussion section. A typical run-time of the code is about 1 min on a standard PC. Here c is the oxygen concentration, j is the proton current density, j e is the electron current density, and η is the local ORR overpotential (positive by convention). Note a small variation of the overpotential η and the oxygen concentration c through the CCL depth and the linear shape of the oxygen concentration through the GDL.

Model
Basic equations.-Schematic of the cathode side of a PEMFC is depicted in Figure 1. The model for the CCL performance has been discussed in detail in Ref. 19 where ε is the Newman's dimensionless reaction penetration depth and μ is the dimensionless constant: The dimensionless variables are defined according tõ are the characteristic time and current density, respectively, and Z is the impedance. In Eqs. 3-5, x is the distance through the CCL depth with x = 0 located at the membrane interface (Figure 1), t is time, j and η are the local proton current density and the ORR overpotential, respectively, b is the ORR Tafel slope, c and c re f are the local and reference oxygen concentrations, respectively, σ p is the CCL proton conductivity, l t is the CCL thickness, C dl is the volumetric double layer capacitance (F cm −3 ), D is the effective oxygen diffusion coefficient in the CCL, and i * is the volumetric exchange current density (A cm −3 ).
Eq. 1 results from substitution of the Ohm's law for the proton current density to the charge conservation equation. Eq. 2 is the oxygen mass balance equation, which implies Fick's law for the oxygen diffusion in the CCL. The right side of Eqs. 1 and 2 is the Tafel law for the oxygen reduction reaction (ORR) rate. Note that the macrohomogeneous model 1, 2 ignores the potential loss due to the oxygen transport in agglomerates of Pt/C particles. Recently, it has been shown that in state-of-the-art catalyst layers with the radius of agglomerate less than 100 nm, this transport is important at high currents only. 27 For the small current densities considered in this work, this transport loss is negligible.
The oxygen transport in the GDL is described by a diffusion equation whereD b is the dimensionless oxygen diffusion coefficient in the GDL, normalized as indicated in Eq. 5. To avoid misunderstanding, the dimensionless oxygen concentration in the GDL is denoted as The system 1, 2, and 7 describes the concerted performance of the CCL and the GDL in a fuel cell. Of key importance are the boundary conditions for this system, as discussed in the next section.
Boundary and interface conditions.-At the CCL/membrane interface, the oxygen flux is zero: The conditions at the CCL/GDL interface (x = 1, Figure 1) express continuity of the oxygen concentration and flux and fix the value of overpotentialη 1 and zero proton current: Generally, the overpotential can be fixed at eitherx = 0, orx = 1; the respective valuesη 0 andη 1 are uniquely related to the cell current densityj 0 . However, fixing the overpotential atx = 1 greatly simplifies solution of the linearized problem for the perturbation amplitudes (see next section). Note that the subscripts 0 and 1 mark the values at the membrane/CCL and the CCL/GDL interface, respectively ( Figure 1). At the GDL/channel interface (x = 1 +l b , Figure 1), the oxygen concentration is equal to the concentration in the channel: [11] wherel b = l b /l t is the dimensionless thickness of the GDL, and c h = c h /c re f is the normalized oxygen concentration in the channel.
Equations for the perturbation amplitude.-Now we can apply small-amplitude harmonic perturbations to the system 1, 2, and 7: where the superscripts 0 and 1 indicate the steady-state solution and the perturbation amplitude in theω-space, respectively. Hereω = ωt * is the dimensionless circular frequency of the exciting signal. Substituting Eqs. 12 into 1, 2, 7 and performing the standard linearization procedure, 19,28 we arrive at a system of linear equations for the perturbation amplitudes ) unless CC License in place (see abstract). ecsdl.org/site/terms_use address. Redistribution subject to ECS terms of use (see 134.94.122.242 Downloaded on 2015-04-28 to IP Solution to the system 13-15 is subjected to the boundary conditions, which follow from linearization of Eqs. 8-11: where we introduced notationÑ 1 1 for the flux of the oxygen concentration perturbation atx = 1. Note thatη 1 1 is the amplitude of the overpotential perturbation applied atx = 1. The system 13-15 is linear, and hence the overpotential can be perturbed at either side of the CCL not changing the results. The impedanceZ of the system is given byZ The system 13-16 describes the "CCL+GDL" impedance at an arbitrary cell current density, provided that the steady-state shapes of the overpotentialη 0 and of the oxygen concentrationc 0 are known. Unfortunately, with the functionsη 0 (x) andc 0 (x) corresponding to medium and large cell currents (see e.g., Ref. 28), analytical solution to the system 13-15 seemingly does not exist. Nonetheless, if the cell current is small, this system has a closed-form solution.
Small current impedance.-If the cell current is small (an exact criterium will be given below), the undisturbed shapes of the oxygen concentration and overpotential are constant alongx: wherec 0 1 is the steady-state oxygen concentration at the CCL/GDL interface. The second of Eqs. 18 is the Tafel law relating the static overpotentialη 0 =η 0 to the cell current densityj 0 (Figure 1).
From a simple balance of fluxes it follows, thatc 0 1 is related to the oxygen concentration in the channelc h viã [19] where [20] is the limiting current density due to the oxygen transport in the GDL. Substituting Eqs. 18 into the system 13-15 and dividing the resulting equations by ε 2 , we come to where =ω ε 2 [24] is the reduced circular frequency. 29 Eqs. 21-23 form a system of linear equations with constant coefficients. However, an attempt to directly solve this system with the boundary conditions 16 using the mathematical software Maple results in an equation of a length of 18 Mb (18 million symbols). A much more compact solution is obtained by using the following strategy. First, we solve Eq. 23, which gives the perturbation amplitude of the oxygen concentration atx = 1 This equation determines the boundary condition for the oxygen concentration perturbation in the catalyst layer, see Eq. 16. Now Eqs. 21 and 22 can be solved; with the solutions forc 1 (x) andc 1 b (x) at hand, the continuity of the oxygen flux atx = 1 allows us to calculate the parameterÑ 1 1 , which completes the solution process. The expressions forc 1 (x) andη 1 (x) are very cumbersome and they are not displayed here. Calculation of the system impedance with Eq. 17 leads to the following result  [26] where the coefficients and arguments are listed in Appendix.
Eq. 26 is valid provided that the cell current density obeys to where j * and j D are the characteristic current densities for the proton and oxygen transport in the CCL, respectively. Physically, Eq. 27 provides constancy of the local overpotential and oxygen concentration through the CCL depth. 26,28,30 It is advisable to calculate the zero frequency limit of the impedance 26. Passing to the limit → 0, expanding the result overj 0 and retaining leading terms only, we get the steady-state differential resistivityR c of the cathode side in the limit of small cell current densityR In the dimensional form this equation reads The first and second terms on the right side of Eq. 29 are the CCL proton resistivity and the faradaic resistivity, respectively. The third and fourth terms are the resistivities due to the oxygen transport in the CCL and the GDL, respectively.

Numerical Results and Discussion
The effect of the oxygen diffusion coefficient in the GDL on impedance spectra is illustrated in Figure 2. The spectra in this Figure are plotted in the coordinatesj 0Z , as suggested in Ref. 10. In these coordinates, in the absence of oxygen transport limitations (D,D b → ∞), the right intercept of the spectrum with the real axis is close to unity. This is immediately seen if we multiply Eq. 28 byj 0 and take into account thatj 0 1. For largeD b = 100, the zero-frequency point of the smaller spectrum in Figure 2a is, indeed, close to unity; the excess is attributed to the sum of the non-faradaic terms in Eq. 28 multiplied byj 0 = 0.1 (Figure 2). The larger spectrum in Figure 2a corresponds to lower D b = 10, and it exhibits two arcs, of which the low-frequency one is due to the oxygen transport in the GDL and the high-frequency arc is due to the faradaic and transport processes in the CCL. The inset in Figure 2a shows the high-frequency part of the spectra, which exhibit the same straight 45 • -line due to the proton transport in the CCL. 16,19 In the standard Nyquist coordinates (Z re ,Z im ), projection  of this straight line onto the real axis gives the dimensionless proton resistivity of the CCLR p = 1/3. In the coordinates of Figure 2, this projection is the proton transport overpotential, which equalsj 0 /3 ( 0.0333 forj 0 = 0.1, Figure 2a). Figure 2b shows what happens to the spectra when the cell current increases. As can be seen, larger current dramatically increases the GDL arc due to much higher oxygen transport losses in the GDL. It should be emphasized that in Eq. 29, the GDL transport resistivitỹ l b /D b is independent of the cell current merely because the respective current-dependent terms have been truncated upon derivation of this particular equation. The resistivityl b /D b is a low-current limit of the generally current-dependent GDL resistivityR G DL , which can be derived from the following arguments.
Omitting the superscript 0, the static overpotential due to the oxygen transport in the GDLη G DL is given bỹ [30] This equation is obtained by substitution of Eq. 19 into Eq. 18 and settingc h = 1. Differentiating Eq. 30 overj 0 , we get a general relation forR G DL :R From Eq. 31 it is seen, thatR G DL l b /D b ifj 0lb /D b 1. However, when the termj 0lb /D b approaches unity, the resistivityR G DL tends to infinity and this trend is seen in Figure 2b. Thus, the contribution  Table I. of the GDL transport to impedanceZ , Eq. 26 is not limited by the conditionj 0lb /D b 1. Eq. 26 is sufficiently "fast" for fitting the experimental Nyquist spectra. By fast we mean any fitting algorithm, which returns the result in less than 10 minutes on a standard PC (a wall-clock time required for a cup of coffee). Figure 3 shows the impedance 26 fitted to the experimental Nyquist spectra of a high-temperature PEM fuel cell. 31 Calculations have been performed in Maple environment; the real and imaginary part ofZ , Eq. 26 have been separated using the Maple commands Re(Z ) and I m(Z ). The fitting was based on the matrix form of the Maple procedure NonlinearFit. Overall, fitting of a single spectrum takes 1-2 min on a standard PC, provided that the initial guess for the fitting parameters is not far from the final values. The Maple worksheet, which implements this algorithm can be downloaded at https:// dl.dropboxusercontent.com/u/79521302/FittingPEMFCspectra.mw The set of fitting parameters for each curve is listed in Table I. The MEA parameters required for fitting are given in the caption to Table I; the cell operating parameters and other relevant experimental information is reported in Ref. 31.
As can be seen, the ORR Tafel slope exhibits 20%-increase as the cell current density changes from 40 mA cm −2 to 120 mA cm −2 ( Table I). The average value of b 50 mV agrees with the number obtained for the same cell from a different impedance model. 22 The CCL proton conductivity also increases with the cell current with the Table I  average of σ p 0.05 −1 cm −1 , which is a typical value for the cell utilizing phosphoric acid as a proton conductor. 32 The double layer capacitance practically does not change with j 0 and the value of C dl 20 F cm −3 correlates with the numbers reported in literature for low-temperature PEM fuel cells. 18 The oxygen diffusion coefficient in the CCL grows linearly with the cell current, with the average of D 2 · 10 −4 cm 2 s −1 . The latter value is close to D = 1.5 · 10 −4 cm 2 s −1 obtained in Ref. 33 using the experimental spectra 21 for the same type of cells. Note that in Ref. 32, D has been determined from analysis of the mass-transport (high-current) region of the cell polarization curve, which yields D 10 −3 cm 2 s −1 . This does not contradict to the numbers in Table I if we extend the linear growth of D with j 0 to the region of large cell currents. This growth can be attributed to increase in the local CCL temperature with the current density.
The last parameter given by this model is the GDL oxygen diffusivity D b . For D b fitting gives a value in the range of 0.02-0.04 cm 2 s −1 (Table I), which roughly is ten times less than the oxygen diffusion coefficient in free air. However, taking into account that in a HT-PEMFC, the GDL is full of phosphoric acid 34 reducing the effective GDL porosity, the numbers in Table I seem to be reasonable. Note that the impedance studies of Bultel et al. 11 gave the numbers in the same range for the oxygen diffusivity in the GDL of a low-temperature PEMFC, which is a similar porous structure partially flooded by water. For the GDL of a low-temperature PEMFC, Malevich et al. 35 reported D b 0.01 cm 2 s −1 , though this number was obtained from the transmission line modeling of the cell impedance. Here we do not discuss the temperature dependence of D b , which is expected to obey to the power law, rather than to the exponential (Arrhenius) law. 36 Table I allows us to verify validity of the model for the largest cell current density of 120 mA cm −2 , which must obey to Eq. 27. Taking for the estimate the oxygen concentration c 1 6 × 10 −6 mol cm −3 (O 2 concentration in air at 150 • C), with the data from Table I we get Thus, the criterium j 0 j * holds for all the cell currents indicated in Figure 3, while for j D we have j 0 j D instead of j 0 j D . This means, that the through-plane shape of the static overpotential corresponding to the semicircle in Figure 3c is not constant, which is not taken into account in the model. Nonetheless, we see that even in this case, the model gives the numbers which are close to those obtained for the lower currents.
The model above ignores the variation of the steady-state shapes of the overpotentialη 0 and the oxygen concentrationc 0 alongx. However, it takes into account the x-variation of the perturbations η 1 andc 1 . Though at small currents, the contributions of the proton conductivity and the oxygen diffusivity to the total impedance 26 are small, the accurate least-squares fitting allows us to obtain the respective transport parameters σ p and D.

Acknowledgment
The author is grateful to Olga Shamardina, Michail Kondratenko and Alexander Chertovich (M.V.Lomonosov Moscow State University) for the experimental impedance data shown in Figure 3.

Appendix: Coefficients and Parameters in Eq. 26
The coefficients and parameters appearing in the expression for the cathode side impedance, Eq. 26 are given by α cc = 2 2y p y mDc α ss = √ 2Dc 0 1j0 −(y m + y p )( p − t) The dimensionless oxygen flux at the CCL/GDL interface is given bỹ where and we have introduced auxiliary parameters Assuming that the oxygen concentration in the channel is constant (high oxygen stoichiometry), the oxygen concentrationc 0 1 is (cf. Eq. 19). The concentrationc 1 1 ≡c 1 b (1) is given by Eq. 25.

Marks dimensionless variables b
Tafel slope (V) C dl Double layer volumetric capacitance (F cm −3 ) c Oxygen molar concentration in the CCL (mol cm −3 ) c 1 Oxygen molar concentration at the CCL/GDL interface (mol cm Oxygen diffusion coefficient in the GDL, cm 2 s −1 F Faraday constant j Local proton current density (A cm −2 ) j 0 Cell current density (A cm −2 ) j * Characteristic current density for proton transport in the catalyst layer (A cm −2 ), Eq. 27 j D Characteristic current density for oxygen transport in the catalyst layer (A cm −2 ), Eq. 27 i Imaginary unit i * Volumetric exchange current density (A cm −3 ) l b Gas-diffusion layer thickness (cm) l t Catalyst layer thickness (cm) Small-amplitude perturbation Greek μ Dimensionless parameter, Eq. 4 ε Newman's dimensionless reaction penetration depth, Eq. 3 η Local ORR overpotential, V η 0 Total ORR overpotential, V η 1 ORR overpotential at the CCL/GDL interface, V σ p CCL ionic conductivity ( −1 cm −1 ) Dimensionless reduced circular frequency, Eq. 24 ω Circular frequency, (s −1 ; ω = 2πf ).