Analytical impedance of oxygen transport in the channel and gas diffusion layer of a PEM fuel cell

Analytical model for impedance of oxygen transport in the gas--diffusion layer (GDL) and cathode channel of a PEM fuel cell is developed. The model is based on transient oxygen mass conservation equations coupled to the proton current conservation equation in the catalyst layer. Analytical formula for the"GDL+channel"impedance is derived assuming that the oxygen and proton transport in the cathode catalyst layer (CCL) are fast. In the Nyquist plot, the resulting impedance consists of two arcs describing oxygen transport in the air channel (low-frequency arc) and in the GDL. The characteristic frequency of GDL arc depends on the CCL thickness: large CCL thickness strongly lowers this frequency. At small CCL thickness, the high-frequency feature on the arc shape forms. This effect is important for identification of peaks in distribution of relaxation times spectra of low--Pt PEMFCs.


I. INTRODUCTION
Electrochemical impedance spectroscopy (EIS) provides invaluable information on transport properties of PEM fuel cell in a current production mode, without interruption of cell functioning 1 . Not surprisingly, EIS of PEM fuel cells is a rapidly growing field 2 . Understanding impedance spectra requires modeling. Strong criticism of equivalent circuit approach has been published by Macdonald in his seminal paper 3 and in recent years, physics-based models for PEMFC impedance tend to replace equivalent circuit modeling (see recent reviews of Tang et al. 2 and Huang et al. 4 ).
Every transport and kinetic process in a PEM fuel cell has its own resonance frequency. If these frequencies do not overlap, one could identify them using the distribution of relaxation times (DRT) technique [5][6][7] . In addition, DRT analysis of impedance spectra returns the contribution of every process into the total differential resistance of the cell. However, correct identification of DRT peaks is a non-trivial task requiring modeling and experimental work. Analytical models predicting characteristic frequencies of oxygen transport processes in the cell could be very helpful in this respect.
Oxygen reduction reaction (ORR) is usually responsible for a large part of potential loss in a PEMFC. Oxygen is transported to the catalyst sites through the cathode channel and gas-diffusion layer; both the transport processes have their signatures in the EIS spectra. After pioneering experimental work of Schneider et al. 8,9 , Kulikovsky and Shamardina 10,11 , Maranzana et al. 12 and Chevalier et al. 13 developed numerical and analytical models incorporating a) ECS Active member; Electronic mail: A.Kulikovsky@fzjuelich.de b) Also at: Lomonosov Moscow State University, Research Computing Center, 119991 Moscow, Russia channel impedance. Formulas for pure channel impedance have been obtained 13,14 assuming fast oxygen transport through the GDL and cathode catalyst layer (CCL). However, the coupling between the channel and GDL impedance remained poorly understood. Recently, Cruz-Manzo and Greenwood reported analytical model for the GDL+channel impedance 15 . However, their result missing important effect of double layer charging on this impedance, as discussed below.
In this work, we develop analytical model for the GDL+channel impedance Z gdlc in a PEM fuel cell. Assuming fast oxygen and proton transport in the CCL, analytical expression for Z gdlc is derived. We show that for typical PEMFC parameters, the Nyquist spectrum of Z gdlc consists of two arcs corresponding to oxygen transport in the channel and GDL. GDL impedance differs from the Warburg finitelength impedance due to "non-Warburg" factor depending on the superficial double layer capacitance of the electrode. For typical PEMFC parameters, the characteristic frequency f gdl of the GDL impedance Z gdl is close to the Warburg finite-length frequency; however, for larger CCL thickness typical for non-Pt cells, the non-Warburg factor strongly lowers f gdl . In the opposite limit of small catalyst layer thickness, a high-frequency feature on the shape of imaginary part of Z gdl vs frequency forms. Analysis shows that the channel impedance also depends on the double layer capacitance meaning that impedance of all oxygen transport medias in the PEMFC "feel" double layer capacitance of the electrode where the oxygen is transported to. The goal of this paper is to clarify the situation with oxygen transport impedance in the GDL and channel after interesting and useful, but incomplete and rather difficult for understanding work of Cruz-Manzo and Greenwood 15 .

II. MODEL
Schematic of the cell with the straight cathode channel is shown in Figure 1. Our main goal here is analysis of impedance of the GDL+channel oxygen transport system. To simplify the model and to separate the GDL+channel impedance from impedance of oxygen transport in the CCL, we will assume that the latter transport is fast. In addition, for simplicity we will assume that the proton transport in the CCL is also fast. The characteristic frequency of proton transport in the CCL is much higher than the other frequencies in the system and hence this assumption does not affect low-and medium-frequency impedances which are of primary interest in this work.

A. Proton charge conservation equation
The proton charge conservation equation reads where C dl is the volumetric double layer capacitance (F cm −3 ), η is the positive by convention ORR overpotential, j is the local proton current density in the CCL, x is coordinate through the cell cathode counted from the membrane, i * is the volumetric exchange current density (A cm −3 ), c is the local oxygen concentration, c in h is the reference oxygen concentration, and b is the ORR Tafel slope.
Since proton and oxygen transport in the CCL are fast, η and c are nearly independent of x. Integrating Eq.(1) over x from 0 to the CCL thickness l t we come to where η 0 is the ORR overpotential at the membrane surface, j 0 is the local cell current density, and c 1 is the oxygen concentration at the CCL/GDL interface.
To simplify calculations, we introduce dimensionless variables is the characteristic time of double layer charging, z is the coordinate along the cathode channel, L is the channel length, D b is the oxygen diffusion coefficient in the GDL, l b is the GDL thickness, ω is the angular frequency of the applied AC signal, and Z is the local impedance.
(2) takes the form Substituting Fourier-transforms into Eq.(5), expanding exponent in Taylor series, neglecting term with the perturbations product, and subtracting the static equation, we get equation relating the small perturbation amplitudesη 1 ,j 1 andc 1 1 in theω-space: Here, the superscripts 0 and 1 mark the static variables and the small perturbation amplitudes, respectively. Local cathode side impedance at a distancez from the channel inlet is given byZ Dividing Eq.(7) byj 1 , we obtain equation forZ: Suppose that the perturbation of oxygen concentration is zero:c 1 1 = 0. Physically, this is equivalent to fast oxygen transport in all transport medias of the cathode side, including channel. Solution to Eq.(9) is theñ which is impedance of a parallel RC-circuit, where the term iω in denominator describes contribution of the double layer capacitance andc 0 1 eη 0 is the inverse faradaic resistivity of the catalyst layer.
The static oxygen concentrationc 0 1 is related to the channel concentrationc 0 h as is the local limiting current density due to oxygen transport in the GDL. Eq. (11) is solution of the static version of equation for oxygen transport in the GDL, see Eq. (18) below. The concentrationc 0 h varies along the channel; numerical model 16 shows that unless the mean cell current density is small, this variation to a good approximation is linear: Note that Eq. (13) is an approximation 1 the exact shape of c 0 h (z) should be calculated using a numerical model 16,17 .
To simplify calculations, we will ignore the factor 1 − j 0 /j lim in Eq. (11), assuming that j lim is large and the variation of static oxygen concentration across the GDL is negligible. With this, Eqs.(7), (9) transform tõ As discussed above, the non-faradaic oxygen transport contributions to local impedance gives the term withc 1 1 . In order to calculatec 1 1 , we need to consecutively solve equations for oxygen transport in the GDL and channel, as discussed in the next section. For further references we need a total faradaic impedancẽ Z f,tot of the cathode, taking into account variation of oxygen concentration along the channel. Suppose that the cell is divided into N → ∞ virtual segments. Local current in each segment flows in the through-plane direction, hence local faradaic impedances are connected in parallel. Thus, Z f,tot is given byZ Electron conductivity of the cell components is assumed to be large and henceη 0 is independent of the coordinatez. Substituting Eq.(13) into Eq.(10), we get the local faradaic impedanceZ f (z). Calculating integral, from Eq.(16) we findZ

B. Oxygen transport in the GDL
Oxygen transport in the gas-diffusion layer is described by the diffusion equation where c b is the oxygen concentration in the GDL and c h is the oxygen concentration in channel. The left boundary condition for Eq. (18) means that the oxygen flux on the GDL side of the CCL/GDL interface equals the local current density in the cell. This condition agrees with the assumption of fast oxygen transport in the CCL. With the dimensionless variables Eq.(3), Eq.(18) takes the form where the dimensionless parameter µ is given by Eq. (19) is linear and hence the equation for the small perturbation amplitudec 1 b is wherec 1 h is the oxygen perturbation in channel (see below). Solution to Eq. (21) is Setting herex = 1, we get the perturbation amplitude of oxygen concentration at the CCL/GDL interfacec 1 b (1) =c 1 where φ and ψ are auxiliary dimensonless parameters

C. Oxygen transport in channel
To a good approximation, oxygen mass transport in the channel can be described by the 1d + 1d plug flow equation: where h is the channel depth. The right side of Eq.(25) is the oxygen diffusive flux in the GDL at the channel/GDL interface representing oxygen "sink" from the channel.
With the dimensionless variables (3), Eq.(25) reads whereJ is the mean current density in the cell ξ is the dimensionless parameter, and λ is the stoichiometry of air flow Eq.(26) is linear and we can immediately write down equation for the perturbation amplitudec 1 h : Differentiating Eq. (22), we find the fluxD b ∂c 1 b /∂x|x =1+l b and Eq.(29) takes the form wherej 1 is a function of coordinatez and ofc 1 h . To find the explicit dependencej 1 (z) we substitute (23) into Eq. (14); solving the resulting equation forj 1 we get D. Solution procedure and total cathode impedance Substitutingj 1 , Eq.(31), into Eq.(30) and solving the resulting equation we find the amplitude of oxygen concentration perturbation along the channel (32) where the independent ofz coefficients A, B and C are given by Setting in Eq.(32)η 1 =Zj 1 and using the result in Eq.(23) we get a formula forc 1 1 , which linearly depends onZ. Substituting thisc 1 1 into Eq.(15), we obtain a linear algebraic equation for local impedanceZ. Solving this equation, we find local impedance, which includes faradaic and oxygen transport in the channel and GDL processes Note that D loc contains complex exponent exp Bz/(λJ) leading to oscillations ofZ along the channel coordinatez (Ref. 11 ). The total cathode side impedance is given bỹ Calculation of integral leads tõ The closed form of integral in Eq.(38) is a key point leading to analytical formula forZ tot . At high frequencies of the AC signal,c 1 h and local impedanceZ are rapidly oscillating functions ofz, which requires a lot of steps in numerical solution of Eq.(29) and makes it difficult numerical calculation of Eq.(38). Analytical result for integral in Eq.(38) solves the problem.

III. RESULTS AND DISCUSSION
Equations of the previous Section contain mean current density in the cellJ and the total static potential loss (overpotential)η 0 . These parameters are related by the polarization curve, which could be obtained from solution of the static version of equations (1), (18) and ( spectra of total cathode side impedance Z tot for the parameters listed in Table I, the mean cell current density J = 0.5 A cm −2 and several stoichiometries of the air flow are shown in Figure 2a. Figure 2b shows the frequency dependence of Im (Z tot ). As can be seen, the impedance consists of three arcs, of which the low-frequency (LF) one strongly depends on λ (Figure 2a). In the pioneering experiments of Schneider et al. this arc has been associated with the oxygen transport in channel 8,9 . The LF arc vanishes as λ → ∞ (Figures 2a,b).
The medium-frequency (MF) arc in Figure 2a represents oxygen transport in the GDL. This arc is not fully seen due to masking effect of the third, high-frequency (HF) faradaic arc (Figure 2). To emphasize the GDL arc, we subtract the total faradaic impedance (17) from the total cathode impedance, Eq.(39). This leads to the GDL+channel impedanceZ gdlc : Nyquist spectra of Eq.(42) for the same parameters exhibit two arcs; now the left, GDL arc is fully resolved (Figure 3a). Frequency dependence of Im (Z gdlc ) is shown in Figure 3b. Two peaks in Figure 3b correspond to the channel and GDL arcs in Figure 3a. At low stochiometry of the air flow, the GDL arc is located at the right "wing" of the channel arc; the latter contributes to imaginary part of the GDL impedance ( Figure 3b). However, as λ increases, the channel arc vanishes and only GDL arc is left (the curves for λ = 1000 in Figures 3a,b). Equation for GDL impedance in the limit of infinite air flow stoichiometry has been derived by Kulikovsky and Shamardina 11 . In the notations of this work, Eq.(36) of Ref. 11 is The frequency-independent factor (1 −Jl b /D b ) in denominator describes the growth of static resistivity upon ap-  11 . In other words, the GDL impedance "feels" the double layer capacitance of the attached CCL, since CCL provides a boundary condition for oxygen transport through the GDL. Similar capacitive correction to the classic Warburg impedance has been discussed by Barbero 20 in the context of Poisson-Nernst-Plank model for the planar electrode. It is worth noting a detailed study of the effect of non-capacitive boundary conditions on Warburg impedance published recently by Huang 21 .
The non-Warburg factor (1 + iω/J) in Eq.(44) leads to important effects. In the dimension form, iω/J = iωC dl l t b/J, i.e., this factor includes the superficial double layer capacitance C dl l t (F cm −2 ) of the CCL. Under constant volumetric DL capacitance C dl (F cm −3 ), this results in dependence of the GDL spectrum on the catalyst layer thickness l t . Figure 4 shows that with typical for Pt/Cbased PEMFCs CCL thickness in the range of 10 µm to 1 µm, the effect of non-Warburg factor on the characteristic frequency f gdl of GDL spectrum is small and f gdl is close to the Warburg finite-length frequency f W : However, for CCL thickness on the order of 100 µm typical for cells with non-Pt catalysts 22 , the non-Warburg factor strongly shifts the frequency f gdl to lower values ( Figure 4). This may lead to partial overlapping of the GDL and channel peaks of the respective imaginary parts (cf. Figure 3b, the curve for λ = 5). Figure 4 also shows that with the decrease in CCL thickness, the width of the GDL peak increases. Moreover, at low l t , the non-Warburg factor leads to formation of a distinct high-frequency feature in the spectrum (Figure 4). The effects in Figure 4 are particularly important for identification of DRT peaks in the high-and low-Pt PEM fuel cell spectra. The thickness of a low-Pt CCL is typically three to four times less than the high-Pt CCL thickness, while the volumetric C dl in both types of cells is the same.
Finally, a formula for "pure" channel impedanceZ chan can be obtained from Eqs.(42), (39) by setting GDL thickness to zero:Z where N c and D c are given in Appendix. The characteristic frequency f chan of channel impedance depends almost linearly of the cell current density and on the air flow stoichiometry ( Figure 5). It is interesting to note that the slope of f chan (J) increases with the growth of stoichiometry λ (Figure 5a), although the magnitude ofZ chan rapidly decreases with λ. due to channel; in the dimension form this resistivity reads R chan = b J 2 (2λ − 1) ln(1 − 1/λ) − λ ln 1 − 1 λ (47) Eq.(47) is independent of parameter ξ, Eq.(28); however, ξ affects the shape of impedance (46). As ξ contains the product C dl l t , the shape ofZ chan appears to be dependent of the CCL thickness. Though for typical PEMFC parameters this dependence is weak, we should stress that impedance of oxygen transport layers (GDL and channel) depend on the double layer capacitance in the electrode.

IV. CONCLUSIONS
An analytical model for impedance of oxygen transport in the GDL and channel of a PEM fuel cell cathode is developed. The model is based on 1d+1d-coupled oxygen mass transport equations along the channel and through the GDL. The oxygen and proton transport through the cathode catalyst layer are assumed to be fast. Analytical solution for the GDL+channel impedance Z gdlc is obtained. The Nyquist spectrum of Z gdlc consists of two arcs, of which the left arcs corresponds to oxygen transport in the GDL and the right (LF) arc is due to oxygen transport in the channel. The characteristic frequency of the channel arc depends linearly on the oxygen stoichiometry. The shape of the GDL arc differs from the Warburg finite-length impedance due to effect of double layer charging in the catalyst layer attached to the GDL. For typical PEMFC parameters, the characteristic frequency of the GDL arc f gdl is close to the Warburg finite-length frequency. However, for larger CCL thickness the non-Warburg factor strongly shifts f gdl to lower values. At low catalyst layer thickness, the high-frequency feature in the spectrum of imaginary part of the GDL impedance forms. This effect is particularly important for identification of DRT peaks in low-Pt cell spectra, as low-Pt cells typically differ from the high-Pt cells by the catalyst layer thickness only.