A Model for Extraction of Spatially Resolved Data from Impedance Spectrum of a PEM Fuel Cell

We report a novel approach to processing of impedance spectra of a PEM fuel cell. We split the cell into N virtual segments and let each segment to have its own set of transport and kinetic parameters. The impedance of a single segment is calculated using our recent physics–based impedance model; the segments are “linked” by equation for the oxygen mass balance in the cathode channel transporting the local phase and amplitude information from one segment to another. Thanks to this transport, the total cell impedance contains information on the local transport and kinetic properties of the cell. We show that ﬁtting the model cell impedance to the experimental spectra yields the parameters of individual segments, i.e., the shape of the cell physical parameters along the cathode the which permits of the work any

Expected wide penetration of polymer electrolyte membrane fuel cell (PEMFC)-based power sources on a market puts forward a problem of fast and reliable characterization of cells and stacks. One of the best in-situ methods for PEMFC characterization is electrochemical impedance spectroscopy (EIS). Impedance spectra contain invaluable information on separate contributions of various processes inside the cell/stack into its performance. However, extraction of this information requires quite sophisticated modeling.
Typical PEMFC operates at an oxygen stoichiometry of two, meaning that the oxygen concentration c at the channel outlet is half of the inlet concentration. Strong non-uniformity of c may induce substantial non-uniformities in the local cell current and, as a result, in the local transport parameters of the cell components. These nonuniformities have been studied using impedance spectroscopy of segmented cells. [1][2][3][4][5][6] However, to the best of our knowledge, the spectra of the non-segmented cell have always been processed assuming uniformity of the cell transport and kinetic parameters over the cell active area.
Of particular interest is the distribution of water over the cell surface. Water is a critically important ingredient for PEM fuel cell operation. Polymer membrane conducts protons only in a wet state and hence without water the cell would "die". On the other hand, it is well known that excess of liquid water produced on the cathode side tends to accumulate in the porous layers, thereby reducing oxygen transport through these layers (flooding). Intelligent water management is, therefore, one of the key issues in the PEMFC technology.
Fuel cell flooding is a complex phenomenon, which includes accumulation of water produced in the oxygen reduction reaction (ORR) in the cathode catalyst layer (CCL) and the gas-diffusion layer (GDL), formation of small liquid water droplets in the cathode channel, merging of these droplets to form a liquid water slug, and transport of the slugs toward the outlet by the pressure gradient. Detailed experimental studies of these steps have been reported 7-10 and a lot of models of various complexity have been developed to simulate the process [11][12][13][14][15][16][17] (see also reviews 18,19 ).
The idea of using impedance spectroscopy for detection of cell flooding is not new. Le Canut et al. reported an experimental study of the effect of PEMFC stack flooding on the impedance spectrum. 20 They showed that flooding increases the impedance at the frequencies below 10 Hz and lowers the impedance phase angle at the frequencies less than 100 Hz. Roy and Orazem suggested to detect the onset of flooding in the cell by statistical analysis of experimental noise of impedance spectra. 21,22 Fouquet et al. have shown that flooding affects the parameters of the Randles-like equivalent circuit fitted to experimental spectra and discussed the respective changes of the circuit parameters. 23 The local oxygen concentration in the channel strongly depends on the local rate of oxygen transport through the porous layers. AC perturbation of the cell potential/current generates local perturbations of the oxygen concentration c 1 everywhere over the cell active surface. However, c 1 is not only produced locally, but also transported along the channel with the flow. Transport leads to mixing of local phases and amplitudes of c 1 with those transported along the channel; in this way, the spectrum of the whole cell contains information about the spatial distribution of the transport coefficients (or any other local parameter) over the cell surface. This gives us a chance to measure the shapes of the cell parameters along the channel using a physics-based impedance model.
Physical segmentation of a cell is difficult to realize; this requires rather sophisticated hardware, multichannel EIS analyzer, etc. Typical fuel cell is not segmented; however, spatial non-uniformities in a cell are of great interest for researchers and developers. In this work, we develop a novel approach to impedance spectra processing. In our model, we split the cell into N virtual segments and we let each segment to have its own set of the transport and kinetic parameters. The impedance model we use for an individual segment has been reported in Refs. 24, 25. According to a quasi-2D approach, the segments are connected by the flow in the cathode channel; the respective equation for the oxygen mass transport along the channel takes into account spatial variation of the cell local parameters. The inverse total cell impedance is then calculated as a sum of inverse impedances of the individual segments. The expression for the total cell impedance can be fitted to the experimental spectrum; in this way, we formulate a least-squares minimization problem with multiple parameters. Below, we report a model with ten and twenty virtual segments and five parameters for each segment, which leads to the least-squares problem with 50 and 100 parameters, respectively. Thanks to extremely powerful minimization algorithms available nowadays, this problem can be solved and the result gives us the spatial distribution of the key cell transport and kinetic parameters along the air channel.

Experimental
The experimental work was performed using a segmented cell system and a test station developed at Hawaii Natural Energy Institute.  The segmented cell approach employs close loop Hall sensors (Honeywell CSNN 191) and data acquisition system which allows us to perform simultaneous measurements of locally resolved current, voltage and impedance. A cell with segmented hardware is operated as a single cell using the test station and standard protocols. Details of the segmented cell system are provided in our previous publications. 24,25 The segmented cell hardware is based on 100 cm 2 cell and contained a standard non-segmented and a segmented flow field plates with the same channel design. The segmented plate consists of ten segments forming a continuous path along ten parallel serpentine channels. Each segment is equipped with its own current collector and GDL and it has an area of 7.6 cm 2 . The segmentation is applicable to either the anode or the cathode; in our tests, the cathode was segmented. The reactant streams were arranged in a co-flow configuration, and segment 1 is the inlet segment, and segment 10 is the outlet.
The experiments were performed with commercially available 100 cm 2 catalyst coated membrane from Gore with Pt/C loading of 0.1 mg Pt cm −2 for the anode and cathode. The thickness of the reinforced membrane was 16-18 μm, while the thickness of the catalyst layers varied in the range of 2-4 μm. Sigracet 25BC (thickness 225-230 μm, 80% porosity) was used as GDLs for both electrodes. 25BC consists of carbon paper substrate and microporous layer with the thickness of 40-45 μm. Segmented GDL was used on the cathode side, whereas a uniform GDL was applied at the anode. The gasket material was made of Teflon, with the thicknesses of 125 μm for the electrodes.
The impedance measurements were performed with H 2 /air conditions and at a cell operating temperature of 80 • C. The anode/cathode conditions were 2/2 stoichiometry, 100/50% relative humidity and 148/148 kPa backpressure. The selected frequency range for the EIS experiments was 0.05 Hz to 10 kHz and the amplitude of the sinusoidal current signal perturbation was 2 A, which resulted in a cell voltage response of 10 mV or lower. Spectra have been measured simultaneously from 10 segments and from the whole cell. For the reasons discussed below, in this work, we used only the spectra from the whole cell for processing. The cell operating and geometrical parameters are listed in Table I.

Impedance Model
Basic equations.-The impedance model used in this work has been discussed in detail in our previous publications. [24][25][26] Here, we briefly list the basic equations and discuss the model extension. Consider a segmented cell shown in Figure 1. The core of the model forms a system of transient equations for the charge and oxygen mass conservation in the CCL of an individual segment: Here, C dl is the double layer volumetric capacitance (per unit electrode volume, F cm −3 ), η is the positive by convention oxygen reduction reaction (ORR) overpotential, t is time, j is the local proton current density, x is the distance from the membrane through the electrode depth, i * is the volumetric exchange current density (A cm −3 ), c is the local oxygen concentration, c in h is its reference (inlet) concentration, b is the ORR Tafel slope, σ 0 is the CCL proton conductivity at the membrane interface, s(x) is the conductivity shaping function, and D ox is the effective oxygen diffusion coefficient in the CCL.
Eq. 1 is the proton charge conservation, Eq. 2 is the Ohm's law relating the proton current density to the gradient of overpotential, and Eq. 3 is the oxygen mass balance equation with the Fick's law of diffusion. Detailed discussion of Eqs. 1-3 is given in Ref. 26.
The CCL proton conductivity is assumed to be exponential function of the coordinate through the CCL depth: where β is the inverse characteristic scale of the exponent, and l t is the CCL thickness. From Eqs. 1-4 it is seen that CCL performance is governed by six parameters: C dl , i * , b, σ 0 , β and D ox . In a cell with uniform over the active area parameters, the exchange current density i * cannot be measured by EIS, as variation of i * shifts the cell polarization curve as a whole along the potential axis, not changing the slope of the curve. In our model, i * has been fixed, although in future calculations we plan to make it variable. To further accelerate the fitting procedure, the parameter β has been determined from analysis of a high-frequency (HF) part of the spectrum at small cell current density, and it has been fixed at β = 5. Physically, the HF part of the Nyquist spectra of a cell with uniform proton conductivity should be a straight line with the 45 • slope. 27 However, spectra of the Gore MEA used in this work exhibit a quasi-linear HF shape with the slope strongly differing from 45 • ; this shape can be explained by non-uniform through the CCL depth proton conductivity. [28][29][30][31] The exponential function (4) provides a good fit of the HF part of the model spectra to experimental impedance data. The remaining four parameters C dl , b, σ 0 , and D ox describing the CCL function are declared below as fitting ones.
The oxygen mass transport in the GDL is also described by the Fick's type diffusion equation where c b is the oxygen concentration and D b is the oxygen diffusion coefficient in the GDL. The segment problems are "linked" by the oxygen mass transport in the channel: where c h is the oxygen concentration in the channel, v is the flow velocity (assumed to be constant), z is the distance along the channel F293 counted from the inlet, h is the channel depth, and is the oxygen flux in the GDL at the CDL/channel interface, where l b is the GDL thickness. Solution to Eq. 6 provides a boundary condition for the local segment problems. On the other hand, the flux N b , Eq. 7, is related to the oxygen flux through the CCL/GDL interface, and thereby it is related to the rates of proton transport and faradaic processes in the individual segments. In this way, Eq. 6 transports the information on the local cell properties down the cathode channel.
The system of Equations 1-7 is completed with the boundary and interface conditions, which express continuity of the oxygen concentration and flux at the channel/GDL and GDL/CCL interfaces, zero proton current at the GDL/CCL interface, and zero oxygen flux at the membrane interface. 24,25 In this work, we have used a more accurate version of the model, 25 which employs the exact expression for the oxygen fluxes at the CCL/GDL interface (see below).
Cell segmentation.-The model Equations 1-6 contain five fitting parameters: b, C dl , σ 0 , D ox and D b . The system of Equations 1-6 has been nondimensionalized, linearized and Fourier-transformed to yield a system of linear equations for the perturbation amplitudes of the overpotential and oxygen concentration in the ω-space. 24,25 Suppose that the cell is separated into N virtual segments and in each segment, we have a unique local set of parameters discussed above. The segments are linked by the flow in the channel, which transports boundary condition for the perturbation of the oxygen concentration to each segment (Figure 1). To calculate the system impedance, we have to solve a system of "local" equations for the perturbation amplitudes in each segment, plus equation for the transport of the oxygen concentration perturbation in the channel. This system is solved numerically in Python environment using the SciPy library. First, a Fourier-transformed version of Eq. 6 is solved, assuming fast oxygen transport in the CCL. Under this assumption, the flux on the right side of this equation has analytical solution, which is independent of the oxygen flux at the CCL/GDL interface. Second, we set N = 10 and in every segment we solve the local through-plane problem along x, which yields the local impedance Z seg,n = − η 1 σ 0,n ∂η 1 /∂ x x=0,z=zn [8] where η 1 is the amplitude of overpotential perturbation, and σ 0,n is the proton conductivity at the CCL/membrane interface in the nth segment. The local solution allows us to correct the flux in Eq. 6 and to repeat the procedure above. These iterations are performed until the desired accuracy of the local impedance is achieved. Finally, since the segment impedances are connected in parallel, the total cell impedance Z cell is calculated as

Results and Discussion
Eq. 9 has been fitted to the experimental spectra using the trust region reflective (TRF) algorithm for nonlinear least-squares fitting procedure from the Python library SciPy. The description of this algorithm states that "The algorithm iteratively solves trust-region subproblems augmented by a special diagonal quadratic term and with trust-region shape determined by the distance from the bounds and the direction of the gradient. This enhancements help to ... efficiently explore the whole space of variables. 32 " Our experience shows that depending on the initial conditions, the TRF method converges to solutions having only minor differences, i.e., the method finds solutions close to the global minimum of the merit function. For the problem The shapes of the oxygen diffusion coefficient in the CCL along the air channel calculated with the fitting parameters for the spectra in (a). Also shown is the shape of D ox calculated with 20 segments. with ten segments we have fifty fitting parameters, five per each segment. To reduce the time of calculations a parallel version of the code has been developed using the MPI library; on a PC with the 2.4 GHz quad-core processor fitting of a single spectrum typically takes about 30 minutes. The advantage of parallel code is that it allows us easily change the number of segments not changing a single line in the code. For comparison, below the curves obtained with 20 segments (100 fitting parameters) are shown.
The experimental and fitted model spectra with 10 virtual segments are shown in Figures 2a-5a. The MEA used in this work was the same, as in Ref. 25; however, the spectra used in Ref. 25 have been measured at large oxygen stoichiometry λ, when the effects of flooding are negligible and the approximation of independent of z cell parameters is expected to work well. Here, the spectra measured at λ = 2 have been taken for processing.
The quality of spectra fitting is good (Figures 2a-5a). Figures 2b-5b show the z-shapes of the oxygen diffusion coefficients in the CCL D ox obtained with 10 and 20 virtual segments. The other parameters do not exhibit strong variation along z; the dependence of mean over the cell surface values of these parameters are shown in Figure 6.
At J = 100 mA cm −2 , D ox is strongly non-uniform: it decays from the channel inlet to outlet by a factor of three (Figure 2b). Note that the shapes of the curves D ox (z) obtained with 10 and 20 virtual segments are close to each other. Both calculations with 10 and 20 segments have started from the uniform distribution of D ox along the channel; close final shapes of D ox (z) indicate that the model works well.
Qualitatively the same behavior demonstrates D ox at the current density of 200 mA cm −2 ; however, with the growth of J , the amplitude of D ox variation decreases and at the cell current of 400 mA cm −2 , D ox is constant along the channel (Figures 2b-4b). Some 20%-amplitude "jumps" of D ox at J = 800 mA cm −2 (Figure 5b) are seemingly  artefacts introduced by the fitting algorithm. Note that the absolute value of D ox increases with J (Figures 2b-5b); a feature that has been reported for the standard high-Pt loaded electrodes. 33,34 The origin of this growth yet is not clear; one may speculate that with the growth of cell current, the temperature and pressure gradients form, which accelerate liquid water removal from the CCL.  The curve for the 20 segments has not been obtained due to poor convergency of numerical procedure.
The strong variation of D ox along z at J = 100 mA cm −2 can be explained by accumulation of water in the CCL close to the channel outlet. The experiments have been performed at the constant oxygen stoichiometry of λ = 2. Under this condition, the flow velocity in the cathode channel increases proportionally to the cell current. We presume that higher flow velocity facilitates liquid water removal from the porous layers, which makes the CCL oxygen diffusivity uniform over the cell surface at higher cell currents (Figures 4b, 5b). At low  currents, the rate of water removal is low and part of the CCL close to the channel outlet appears to be flooded (Figure 2b).
In contrast to D ox , the parameters b, C dl , σ 0 and D b are nearly constant along the channel. The first three parameters in this list (b, C dl and σ 0 ) depend on nanoscopic details of Pt/electrolyte interface in the electrode. Their independence of electrode flooding means that flooding only happens in large pores of the CCL, through which gaseous oxygen is transported. Excess of liquid water in large pores does not change the situation in small pores inside the Pt/C agglomerates, which is responsible for the properties of Pt/Nafion interface. In other words, Pt/C agglomerates remain fully "flooded" regardless of the amount of liquid water in large pores.
On the other hand, b and σ 0 exhibit linear growth with the cell current density, as Figure 6 shows. The growth of b is seemingly due to "cleaning" of Pt surface from oxides with lowering of the cell potential. 35 The growth of the proton conductivity with the cell current shown in Figure 6b has been detected also in high-Pt electrodes. 36 This feature could be explained by the growth of a number of protonconducting channels in Nafion cluster at higher current densities.
The model discussed provides a useful tool for analysis of PEMFC impedance spectra. Thanks to the transport of amplitude and phase information with the flow in the channel, the spectra of the whole cell contain information on spatial distributions of the key parameters over the cell surface. The model above allowed us to extract these distributions. Additional arguments in favor of this idea based on equivalent transmission line approach are discussed in Appendix.
In 2007, Schneider et al. published pioneering works devoted to experimental study of the channel transport impedance in PEMFCs. 2,3 They have shown that this impedance manifests itself as a separate low-frequency arc in the impedance spectra and clarified the effect of oxygen transport in the channel on local segment spectra. Here, we show that the role of the channel impedance is even more important: it allows us to extract the spatial distribution of parameters over the cell surface.
In this work, the spectra from individual segments have been ignored and only the spectra from the whole cell have been taken for processing. "Local" spectra have been measured applying a synchronous AC signal to all the segments. This means that every local spectrum contains information on the transport properties of the upstream domain. Simultaneous fitting of the local spectra could, in principle, be done using the similar approach. This, however, poses significant programming problems; this work is in progress.

Conclusions
We report "segmented" physical model for PEMFC impedance, which includes the impedance due to oxygen transport in the channel. The cell is separated into N virtual segments; a key feature of this model is that each segment has its own set of the transport and kinetic parameters. Thanks to the air flow in the channel, the information on local segment properties is transported down the channel by the phase and amplitude of the oxygen concentration perturbation. This means that the total cell impedance contains information on the spatial distribution of the cell parameters along the channel. The total impedance of this system is calculated using our recent physics-based model. 24,25 Fitting the model impedance to the experimental spectra yields the spatial variation of the cell parameters along the channel. At the cell current densities of 100 and 200 mA cm −2 , the oxygen diffusion coefficient in the cathode catalyst layer exhibits strong decay toward the channel outlet. This behavior could be attributed to the insufficient rate of liquid water removal from the cell.  thank Gunter Randolf for valuable support regarding the system operation.

Appendix: Transmission Line Analogue
To elucidate the role of oxygen transport in the channel, consider first the transmission line (TL) depicted in Figure A1. This system is a crude model of a segmented fuel cell with ideal transport of oxygen in the channel. For simplicity, the impedance of each segment is represented by an RC-circuit. Suppose that due to some spatial non-uniformities, the resistivities and capacitances of the segments are different. Is there a chance to locate the position of each RC-circuit in this line by means of impedance spectroscopy? Evidently, the answer is no; the total inverse impedance 1/Z N of the system with N elements is a sum of inverse RC-impedances and permutation of the terms in this sum does not change the result. Consider now the TL in Figure A2a, where r n model the resistivities due to oxygen transport in the channel. To calculate the impedance of this system, we redraw it in the form shown in Figure A2b. From this Figure it follows that the impedance of the system with N elements is given by recurrent relation where Z N −1 is the impedance of the remaining system with N − 1 elements, and Z RC,1 is the impedance of RC-circuit given by the standard expression Replacing N → N − 1 in Eq. A2 and substituting the resulting Z N −1 into Eq. A2, we get Setting now N → N − 2 in Eq. A2 and using the result in Eq. A4, we come to This process of substitutions can be continued; we see that the total impedance of the system in Figure A2a is given by continued fraction of the form Eq. A5. In this expression, the impedance of each RC-circuit (segment impedance) occupies a separate unique level in the right-side fraction. Clearly, if we permute two RC-circuits in Figure A2a, the respective impedances will exchange their levels in Eq. A5, which will ) unless CC License in place (see abstract). ecsdl.org/site/terms_use address. Redistribution subject to ECS terms of use (see 207.241.231.83 Downloaded on 2019-04-28 to IP inevitably change the total impedance Z N . Thus, the presence of the resistivities r n in Figure A2a makes the total system impedance uniquely dependent on the position of each RC-element in the system.