A Fast Low-Current Model for Impedance of a PEM Fuel Cell Cathode at Low Air Stoichiometry

We report a physics-based analytical model for low-current PEM fuel cell impedance, which takes into account oxygen transport in the channel. The model is fast: ﬁtting the model impedance to the experimental Nyquist spectrum takes about 2 minutes on a 2-GHz notebook. In addition to the transport and kinetic parameters of the catalyst layer, ﬁtting returns the oxygen diffusion coefﬁcient in the gas-diffusion layer and the resistivity due to oxygen transport in the channel. Maple worksheet with the ﬁtting procedure is available for download.

One of the best in situ methods for fuel cell characterization is electrochemical impedance spectroscopy (EIS). Application of a small AC potential perturbation to a working fuel cell allows one to separate contributions of various transport and kinetic processes running in the cell thanks to difference in time constants associated with the physical processes 1,2 However, understanding experimental impedance spectra requires quite sophisticated modeling. Over recent years, a lot of efforts has been done to develop physical models for a PEMFC impedance. [3][4][5][6][7][8][9][10][11][12][13][14] Most of these models are numerical and they are slow to be used in algorithms for spectra fitting. A "faster" alternative is analytical modeling of impedance; [15][16][17][18][19][20][21][22][23] however, fast analytical models have been developed for the case of infinite stoichiometry of the air (oxygen) flow in the cathode channel.
A typical PEMFC operating stoichiometry is around 2 and a typical impedance spectrum exhibits a low-frequency arc due to oxygen transport in the channel. 20,[24][25][26] "Cutting" of this arc allows one to fit the remaining points using a 1D through-plane model; however, discarding the low-frequency (LF) points does not allow us to determine the oxygen diffusion coefficient in the gas-diffusion layer (GDL) and the resistivity due to oxygen transport in the channel. 21 In this work, we develop a most general physics-based analytical model for PEM fuel cell cathode impedance at low current density. The model takes into account the impedance due to oxygen transport in the air channel and hence it is applicable for fitting the spectra measured at low oxygen stoichiometries. The model is fast: fitting returns a set of key PEMFC transport parameters in just two minutes on a 2-GHz notebook.
Analytical expressions for the components of the cathode side impedance have been derived in previous works; 20,22,27 however, plain summation of the respective components leads to a very slow fitting code. Here, we derive a "fast" analytical expression for the total impedance of the cathode side and demonstrate this expression "at work" by fitting it to the experimental spectra of a low-and high-temperature PEMFCs. Maple worksheet with the fitting procedure is available for download at https://github.com/akulikovsky/ Fitting_Procs/issues/2

Model
The model includes three impedances: r The faradaic and proton transport impedance, taking into account a non-uniform distribution of proton conductivity through the cathode catalyst layer (CCL) depth. r The impedance due to oxygen transport through the CCL. * Electrochemical Society Member. z E-mail: A.Kulikovsky@fz-juelich.de r The impedance due to oxygen transport in the GDL and channel.
Analytical expressions for all these impedances are derived assuming that the mean cell current density J is small. This means that the following relation must hold Here,σ p is the mean through the CCL depth proton conductivity, b is the ORR Tafel slope, l t is the CCL thickness, D ox is the effective oxygen diffusion coefficient in the CCL, and c 1 is the oxygen concentration at the CCL/GDL interface. Eq. 1 means that J must be much less than the characteristic current densities for the proton and oxygen transport in the CCL. 28 For typical PEMFC parameters (Table  I) and the CCL thickness of 10 −3 cm, Eq. 1 holds for the cell current densities below 100 mA cm −2 .
Faradaic and proton transport impedance.-The expression for the combined charge transfer (faradaic) and proton transport impedance Z ct+ p is: 29 Here, J 0,1 and Y 0,1 are the Bessel functions of the first and second kind, respectively, j 0 is the local cell current density (see below), C dl is the double layer capacitance, σ 0 is the CCL proton conductivity at the CCL/memebrane interface, and ω is the angular frequency of the applied signal (ω = 2πf , f is the regular frequency in Hz). Eq. 2 has been derived assuming that the proton conductivity σ p decays exponentially through the CCL depth:  where β is a dimensionless inverse characteristic scale of the exponent, x is the distance from the membrane.

Impedance due to oxygen transport in the catalyst layer.-
The low-current impedance due to oxygen transport in the CCL Z ox is given by 27 where are the characteristic frequencies, is the dimensionless Warburg-like impedance, and j ox is given in Eq. 1. Eq. 5 is obtained assuming that the undisturbed oxygen concentration and ORR overpotential are nearly constant through the CCL depth. 27 Both assumptions hold when the cell current density obeys to Eq. 1.

Impedance due to oxygen transport in the GDL and channel.-
In this section, the sign tilde marks dimensionless variables and parameters, defined as [8] and the dimensionless frequency is given by Here, l b is the GDL thickness, D b is the effective oxygen diffusion coefficient in the GDL, η 1 1 is the amplitude of the applied potential perturbation, and c in h is the inlet oxygen concentration. The combined local impedance due to oxygen transport in the channel and GDL Z gdl+c is given by the following chain of expressions. 20 Here and the static oxygen concentration at the CCL/GDL interface c 0 1 is where c 0 h is the static oxygen concentration in the channel. The perturbation amplitude of the oxygen concentration at the CCL/GDL interface c 1 1 is given by Eq. 13 is equivalent to Eq. A3 in Appendix in Ref. 20.
The parameter c 1 1 depends on the perturbation of the local oxygen concentration in the channel c 1 h , which is Here,J are the dimensionless mean cell current density and distance along the channel, respectively, L is the channel length, λ is the air (oxygen) flow stoichiometry v in is the inlet air flow velocity, h is the channel depth, and Parameter χ is The local cell current density j 0 and the static oxygen concentration in the channel c 0 h depend on the distance along the channelz: 30 This completely defines Z gdl+c , Eq. 10. From Eqs. 15 and 13 it follows that c 1 h ∼η 1 1 and hence c 1 1 ∼η 1 1 ; thus, the impedance (10) is independent of the applied potential perturbation η 1 1 , provided that η 1 1 is small.
Total impedance of the cathode side.-Eqs. 2, 5 and 10 give the components of the local impedance of the cathode side Z loc : Z loc = Z ct+ p + Z ox + Z gdl+c [22] Note that all the three terms in this sum depend on the distance along the channel through this dependence of j 0 , Eq. 21; in addition, Z gdl+c depends onz also through c 1 h , Eq. 15. Suppose that the cell is segmented into N segments, each one having the local impedance Z loc . Noting that all the segments are connected in parallel and taking the limit N → ∞, for the total impedance of the cathode side we obtain Generally, the integral in Eq. 23 has to be calculated numerically, which dramatically slows down the fitting procedure (see below). However, numerical tests show that a good approximation of this integral can be obtained as following. Consider the equation j 0 = J , or, taking into account 21 Solving this equation forz we get a pointz * , where the local current density equals the mean current density: where f λ is given by Eq. 18. Calculations show that for λ ≥ 2, an accurate approximation of integral 23 is From Eq. 23 we, thus, get Eq. 25 leads to much faster fitting code, as discussed below.

Results and Discussion
The model impedances have been fitted to an experimental spectrum of a PEMFC 21 (Figure 1). Fitting has been performed in Maple environment using a built-in fitting procedure NonlinearFit. To start the fitting process, initial values for the fitting parameters have been specified, which give the model spectrum reasonably close to the experimental one. The merit function was the sum of squares (Re(Z ) − Re(Z ex p )) 2 + (Im(Z ) − Im(Z ex p )) 2 , where the subscript ex p denotes the experimental values. Note that the sum of squares above is minimized for every frequency point ω, which means that optimization is performed simultaneously in the Nyquist and Bode coordinates.
The spectrum has been measured at the cell current density of 100 mA cm −2 , the experimental conditions are indicated in the caption to Figure 1; more details are given in Ref. 21. Note that in Ref. 21, the spectra have been measured from individual segments and from the whole cell; here we use the data for the whole cell. The LF arc in Figure 1 represents the combined oxygen transport in the GDL and cathode channel. Two variants of the model above have been fitted: approximate Eq. 25 (open circles in Figure 1), and Eq. 23 using numerical calculation of integral (dashed line in Figure 1). As can be seen, for both the model variants the quality of fitting is very good; only the region between the LF and faradaic arcs is fitted not quite well (Figure 1). Note that the fitting points of the model Eq. 25 (open circles) are shown for the same frequencies as the measured points (filled circles). It is also seen that approximate Eq. 25 gives the curve which is very close to the "exact" curve from Eq. 23; however, the approximate model is nearly two orders of magnitude faster, than the exact model (2-3 min vs 3 hours on a standard PC). Figure 1b shows the components of the total cathode side impedance. The GDL Z gdl and channel impedance Z c have been separated by setting λ = 1000 in Eq. 10; this gives the "pure" GDL impedance Z gdl . The channel impedance is then obtained as Z c = Z gdl+c − Z gdl . As can be seen, the GDL and channel impedances are quite large. The fitting parameters resulted from the two model variants are listed in Table I. First we note, that both the models give almost the same parameters. Further, the Tafel slope and the mean proton conductivity agree well with the literature data. 31,32 The double layer capacitance is 50% higher, than the value reported for this MEA in Ref. 21. Here, however, we used a more accurate equation for the impedance Z ct+ p , which takes into account a rapid decay of the proton conductivity through the CCL depth, while in Ref. 21 this parameter was assumed to be independent of the distance from membrane. This leads to somewhat higherσ p and C dl , as compared to Ref. 21. The oxygen diffusion coefficients in the CCL and GDL by the order of spectrum Z tot and its components: charge-transfer plus proton transport, Z ct+ p , oxygen transport in the CCL, Z ox , oxygen transport in the GDL, Z gdl and in the channel, Z c . To represent the details, individual contributions of impedance are shifted using the values of the respective static resistivities, e.g., the CCL oxygen-transport impedance is shifted to the right along the real axis by the value of R ct+ p etc. (c) Bode plot of the same data as in (a). magnitude agree with those reported, 21 taking into account that D ox increases with the cell current density. 33 Note that the air stoichiometry λ has been declared as a fitting parameter. The reason is that Eqs. 21 and 15 have been derived for the cell with the straight channel. In experiment, however, the cell with the meander-like channel has been used. 21 In this flow field, oxygen is transported under the rib between two adjacent turns of meander, which homogenizes the distribution of oxygen concentration along the cell surface. Experimental shape of local current along the channel of the meander-like flow field shows formation of plateau of j 0 in the middle of the channel (Figure 4 of Ref. 21). This plateau arises seemingly due to under-rib convection of oxygen and it changes the value of integral in Eq. 23. In our 1D+1D model, this "homogenization" of local current is effectively accounted for by higher air flow stoichiometry.
Another example of fitting is shown in Figure 2, which displays experimental and fitted Eq. 25 spectra of a high-temperature PEM fuel cell. 17 The spectrum has been measured at the cell current density of 60 mA cm −2 , cell temperature of 160 • C, pressure of 1 bar and air flow stoichiometry λ = 2; the other experimental details and MEA properties are given in Ref. 17. As can be seen, the quality of fitting is very high: the fitted and experimental points shown for the same frequencies practically overlap. The fitting parameters are collected in Table  II. The parameters agree well with those reported in Ref. 19 for this type of cell, except D b , which is three times lower than D b reported in Ref. 19. However, a feature of HT-PEMFC is a large and uncontrollable amount of liquid phosphoric acid in the GDL, which may strongly affect D b in different conditions. Note a significantly larger impedance due to oxygen transport in the CCL in Figure 2b, as com- pared to that impedance in a low-T PEMFC (Figure 1b). This is explained by the five times thicker CCL in the HT-PEMFC.
To conclude we note that the shape of the low-frequency arc in the spectrum strongly depends on the stoichiometry of the air flow, and on the oxygen diffusion coefficient in the GDL. Thus, fitting the spectra using the model above allows us to reliably determine the GDL oxygen diffusivity and the resistivity due to oxygen transport in the channel.

Conclusions
A fast physics-based model for impedance of a PEM fuel cell operated under low stoichiometry of the air flow is developed. The model includes the charge-transfer and the proton transport impedance, and the oxygen transport impedances in the catalyst layer, gas diffusion layer, and in the cathode channel. The model is fast: fitting the model impedance to a single experimental spectrum takes about 2 min on a standard PC. The fitting returns seven parameters: the ORR Tafel slope, the mean proton conductivity of the CCL, the double layer capacitance, the oxygen diffusion coefficients in the CCL and GDL, the effective air flow stoichiometry, and the parameter describing the rate of proton conductivity decay with the distance from the membrane. The effective air flow stoichiometry is introduced to take into account the under-rib oxygen flow in the meander-like flow field.

Acknowledgments
The author is grateful to Dr. Tatyana Reshetenko (University of Hawaii) and to Dr. Mikhail Kondratenko (Moscow State University) for the experimental data and many useful discussions. Small-amplitude perturbation