Polarization Curve of a Non-Uniformly Aged PEM Fuel Cell

: We develop a semi-analytical model for polarization curve of a polymer electrolyte membrane (PEM) fuel cell with distributed (aged) along the oxygen channel MEA transport and kinetic parameters of the membrane–electrode assembly (MEA). We show that the curve corresponding to varying along the channel parameter, in general, does not reduce to the curve for a certain constant value of this parameter. A possibility to determine the shape of the deteriorated MEA parameter along the oxygen channel by ﬁtting the model equation to the cell polarization data is demonstrated.


Introduction
Aging is one of the key problems in polymer electrolyte membrane fuel cell (PEMFC) technology. In particular, stack aging is one of the largest hurdles on the way toward commercialization of a fuel cell car, where a PEMFC stack must stand long-term operation in a harsh environment under variable load.
Over the past two decades, huge efforts have been directed toward studying the aging phenomena in PEMFCs (we refer reader to the review [1]). The cell aging could be caused by numerous chemical and electrochemical phenomena in cell components (bipolar plates, gas-diffusion layers, catalyst layers and membrane). Among the most detrimental are the aging processes in the cathode catalyst layer (CCL), which converts ionic current into electronic current with the aid of oxygen molecules.
A typical CCL aging scenario includes the corrosion of carbon support leading to the detachment of Pt particles from supporting carbon "balls", with subsequent Pt particles agglomeration [2][3][4][5]. This process is seemingly accompanied by change of the Nafion/carbon surface structure [6] and by decreasing of the CCL porosity, which leads to lowering of the CCL proton conductivity and oxygen diffusivity. Note that carbon corrosion accelerates when the cell operates close to the open-circuit potential. A very fast corrosion occurs during the start-stop cycle, when part of the anode channel is filled with hydrogen, and the other part contains air [2].
Fuel cell typically operates at oxygen stoichiometry λ between 1.5 and two. Under these conditions, the distribution of local current density j 0 along the oxygen channel is non-uniform, with j 0 being maximal at the channel inlet [7]. We therefore may assume that the CCL aging runs faster close to the channel inlet; experimental evidence supporting this assumption has been provided in [8].
The polarization curve is a portrait of a fuel cell; measuring this curve is a simple and straightforward procedure. Numerous fuel cell diagnostic methods developed over the past years are much more time-consuming and expensive; some of the diagnostic techniques require cell disassembling. Can we understand by analyzing the cell polarization curve, what aging processes are running in the cell and where on the cell surface these processes peak? A positive answer to this question would greatly simplify aging studies and allow for fast control of working cells and stacks.
In this work, we develop a semi-analytical quasi-2D model of a PEMFC polarization curve. The model is based on an approximate analytical solution to the through-plane conservation equations in the catalyst layer and gas diffusion layer (GDL), linked to the mass conservation equation for oxygen transport in the channel.
Using this model, we analyze the effect of non-uniform (along the channel) aging of the MEA parameters. We show, that the non-uniform aging can be distinguished from the uniform aging of the respective parameter by the inspection of polarization curves. We demonstrate the fundamental possibility to obtain the shape of the "aged" parameter along the channel by fitting the model to the cell polarization data.

Local Polarization Curve
Unless stated otherwise, the equations below are written using the following dimensionless variables: Here, z is the distance along the channel of a length; L ( Figure 1); j is the current density; c is the local oxygen concentration; c 0 h is the reference (inlet) oxygen concentration; η is the overpotential; b is the Tafel slope of the oxygen reduction reaction (ORR) and D is the oxygen diffusion coefficient in the CCL. The characteristic scaling values for current density j * and oxygen diffusivity D * are given by: where σ t is the CCL proton conductivity and l t is the CCL thickness. A local polarization curve of the CCL has recently been derived in [9]. This curve is given by: Here,η 0 is the total potential loss (overpotential) in the CCL,j 0 is the local cell current density, c 1 is the oxygen concentration at the CCL/GDL interface and ε is the Newman's dimensionless reaction penetration depth: where i * is the volumetric exchange current density (per unit CCL volume, A cm −3 ). The first term in Equation (3) takes into account the ORR activation overpotential and the potential loss due to proton transport in the CCL [9]. Note that at large cell currents, this term describes Tafel slope doubling due to poor proton transport in the active layer. The second term in Equation (3) is a potential loss due to oxygen transport in the CCL.
The parameter, β, which appears in Equation (3), is a solution to equation: A reasonably accurate approximate analytical solution to this equation is [10]: A more elegant solution to Equation (5) can be obtained by analytical iterations [9]. This solution, however, leads to much slower code (see below).
A simple linear diffusion equation gives the following relation between the oxygen concentration at the CCL/GDL interface,c 1 , and the local oxygen concentration in the channel,c h : Here: is the limiting current density due to oxygen transport in the GDL and D b is the oxygen diffusion coefficient in the GDL of a thickness, l b . The superscripts, 0 and 1, mark the values at the channel inlet and outlet, respectively ( Figure 1). Substituting Equation (7) into Equation (3), we get a polarization curve of the cathode side, which takes into account oxygen transport in the GDL: Taking into account the definition of dimensionless variables Equation (1), we see that Equation (9) includes the following basic kinetic and transport parameters: the Tafel slope of the ORR, b, the exchange current density, i * , the CCL oxygen diffusivity, D, and proton conductivity, σ t , and the GDL oxygen diffusivity, D b . The model, therefore, allows us to study aging processes, which affect the five aforementioned parameters.

Quasi-2D Polarization Curve
To take into account the effects of non-uniform cell aging, the parameters,D, ε andj 0 lim , in Equation (9) should be considered as functions of the coordinate,z, along the channel. Parametric non-uniformities enhance a non-uniformity of the local current distribution alongz, and thus,j 0 =j 0 (z). In general, oxygen stoichiometry λ could be finite, which means thatc h also varies along the channel. This variation obeys the oxygen mass conservation equation: where J is the mean current density in the cell. Note that the following relation holds: In view of this relation, it is convenient to introduce a stream function: The stream function satisfies the boundary conditions φ(0) = 0 and φ(1) =J; the second of which results from Equation (11).
Integrating Equation (10) overz from zero toz, we can express a formal solution to this equation in terms of φc Substituting Equations (12) and (13) into Equation (9), we come to: where the subscript,z, indicates differentiation overz. A simpler equation, which corresponds to the infinite oxygen stoichiometry, is obtained from Equation (14) by setting λ → ∞: Equation (15) is a first-order ordinary differential equation (ODE) with respect to φ. However, this equation cannot be resolved with respect to φz, which makes a numerical solution of this equation a non-trivial task.
A simpler way to find φ is to differentiate Equation (15) overz, taking into account that ∂η 0 /∂z = 0 (the cell electrodes are equipotential), and that φ, ε,D,j 0 lim and β are functions ofz. This leads to a second-order ODE for φ, which can be directly solved with the boundary conditions φ(0) = 0 and φ(1) =J. Now, the advantage of the model formulation in terms of φ is evident. Indeed, Equation (15) contains bothη 0 andJ; however, these functions are related by the polarization curve, which a priori is unknown. Thus, the solution of Equation (15) would require iterations. In contrast, the equation, ∂η 0 /∂z = 0, withη 0 given by Equation (15), contains onlyJ, and it can be directly solved.
Unfortunately, differentiation of Equation (15) overz leads to a very cumbersome equation, which is not displayed here. However, this equation can easily be generated from Equation (15) by mathematical software (e.g., Maple , Waterloo Maple Inc., Waterloo, ON, Canada). Below, we report numerical solutions to equation: whereη 0 is given by Equation (15). Note that all the results below are obtained assuming infinite oxygen stoichiometry (λ = ∞), a condition that can always be fulfilled in measuring the polarization curve of a degraded cell.
Equation (16) is a second-order ODE with respect to φ(z), which is solved numerically using Maple . With the function φ in hand, we calculate the local current density,j 0 (z), from Equation (12). The overpotential is then calculated by settingz = 0 in Equation (15). Finally, the cell potential is calculated as V cell = V oc − bη 0 , where V oc is the open-circuit potential.

Polarization Curves
Carbon corrosion in the CCL leads to agglomeration of Pt particles that lost contact with the carbon support; corrosion lowers the CCL porosity and reduces the connectivity of Nafion cluster, thereby lowering the CCL proton conductivity. In terms of the model above, carbon corrosion changes the parameters, ε (through the lowering of i * and σ t ) and D. In addition, experiments show that GDL transport properties also worsen with time, which is equivalent to the reduction of the GDL oxygen diffusivity, D b . To understand the effect of non-uniform aging of these parameters on the polarization curve, we will consider their variations separately. The basic set of parameters corresponding to a pristine cell is listed in Table 1.
In a well-designed cell, local current density is typically higher at the oxygen channel inlet; thus, we may expect that the fastest D and i * degradation occurs nearz = 0. On the other hand, liquid water tends to accumulate in the GDL close to the oxygen channel outlet, and thus, the effect of GDL aging is more pronounced close to the outlet. To mimic these three types of cell degradation, we use the following model functions: where D 0 , i 0 * and D 0 b are the base-case (pristine) values. The factor four in the exponent in Equations (17) Figure 2. Qualitatively, aging of the CCL oxygen diffusivity, D, forces the polarization curve to decay faster and makes the decay more smooth (Figure 2a). Non-uniform aging of the exchange current, i * , shifts the curve down along the potential axis and increases the slope of the curve (Figure 2b). The non-uniform decrease of the GDL oxygen diffusivity, D b , manifests itself as the lowering of the limiting current density (Figure 2c). Figure 2 also shows the curves, which are obtained for the case of the uniform aging of the respective parameter. Polarization curves with the uniform and non-uniform aging of D b are practically indistinguishable (Figure 2c). This means that non-uniform aging of the GDL oxygen diffusivity can hardly be detected by the analysis of a cell polarization curve. It can be shown that in the case of ideal proton conductivity and fast oxygen diffusion in the CCL, the limiting current density is simply an average of the local limiting currents over the channel length [11]. This means that the polarization curve of a cell with variable j lim alongz is always equivalent to a polarization curve of such a cell with certain uniform j lim . Figure 2c suggests that this is a more general result, which holds regardless of the transport limitations in the CCL.
However, the curves corresponding to non-uniform aging of D and i * differ from the respective "uniformly aged" curves ( Figure 2a,b). Physically, at low D, the oxygen reduction reaction runs in a small domain close to the CCL/GDL interface, where more oxygen is available [12]. Due to poor oxygen transport, the rest of the CCL experiences "oxygen starvation", and it does not contribute to current conversion. This regime leads to the doubling of the apparent Tafel slope, which dramatically increases the overpotential [12].
A characteristic "indicator" of this regime is the dimensionless reaction penetration depth for the oxygen transport in the CCL given by [12]: Essentially, in the CCL with high oxygen diffusivity, the parameter,l D , is much greater than unity, indicating that all parts of the CCL are easily accessible for oxygen. In the opposite limit ofl D 1, the electrode works in the oxygen-deficient regime, which leads to the doubling of the apparent Tafel slope.
With the data from Table 1, taking c 1 = c 0 h and the cell current density of j 0 = 1 A cm −2 , we getl D 15. This shows that the pristine CCL is highly "transparent" to oxygen. However, with the 55-times lower D in the aged electrode, we getl D 1, which means, that at the channel inlet, the electrode experiences strong oxygen starvation. Further, non-uniform aging of CCL diffusivity leads to the redistribution of local current density alongz, as compared to the uniformly aged cell. Part of the total current is produced in the oxygen-starving domain with twice the Tafel slope, while the other part is produced in the domain with "normal" Tafel kinetics. Figure 2b shows that this situation cannot be described by any constant alongz oxygen diffusivity. The curve for a constantD = 0.1D 0 does not fit well the non-uniformly aged curve (Figure 2b).
Note that much lower oxygen diffusivity in the GDL, D b , could also lower the parameter,l D , by lowering the oxygen concentration, c 1 , at the CCL/GDL interface. However, any variation in D b is easily recognizable in the cell polarization curve as a variation of the limiting current density.
Qualitatively, the effect of non-uniformly aged i * is similar to the effect of non-uniform D. In a cell with a non-uniformly degraded active catalyst surface, part of the current is produced in the low-i * region, where it "costs" more polarization potential. With the growth of the mean cell current, the contribution of the aged domain increases, which leads to the faster decay of the cell polarization curve, is compared to the uniformly aged i * (Figure 2b).
It is worth mentioning that the shapes of the local current alongz corresponding to the non-uniformly aged D and i * are close to each other (Figure 3). This means that measurement of the local cell current (using, e.g., the segmented cell technique) does not guarantee reliable detection of the aging mechanism. ( ) Doubling of the apparent Tafel slope could also be caused by lowering of the CCL proton conductivity. Though this effect is not considered here, it can be understood using the following arguments. The dimensionless reaction penetration depth, l σ , for proton transport in the CCL is given by [12]: The inequalityl σ 1 indicates the regime with the apparent Tafel slope doubling. Physically, in that case, the reaction runs in a narrow domain close to the membrane/CCL interface, where protons are "cheaper". The rest of the CCL is inactive, due to poor proton transport.
With the data from Table 1 and j 0 = 1 A cm −2 , we getl σ 0.9. This means, that the pristine cell works in the intermediate regime between the normal-and double-Tafel regimes. However, with several times lower σ t , the value ofl σ would be much less than unity, which means a double-Tafel regime with much larger polarization loss. Thus, the non-uniform aging of the CCL conductivity cannot, in general, be represented by any equivalent cell with uniform σ t .
The results above are obtained assuming that the oxygen stoichiometry is infinite; this assumption can always be fulfilled in experiments with aged cells. Calculations show that finite λ does not emphasize the difference between uniformly-and non-uniformly aged curves. However, finite λ makes the numerical solution of the problem much more time-consuming.
To summarize, non-uniform aging of parameters D, i * and σ t could, in general, be understood by the analysis of cell polarization curves. Such an analysis requires accurate fitting of the "aged" polarization curves using the quasi-2D cell model, as discussed in the next section.

Can We Understand Non-Uniform Aging by Polarization Curve Fitting?
To address this question, the following procedure has been performed. We generated sixteen points belonging to the non-uniformly aged polarization curve in Figure 2a (exponential alongz aging of the CCL oxygen diffusivity). The cell potential was then perturbed by adding a random δV in the range −10 mV ≤ δV ≤ 10 mV, which mimics experimental errors when measuring the polarization curve.
We assumed that the pristine value ofD 0 is known, e.g., from fitting the pristine cell polarization curve, and we sought the "unknown"z-shape ofD(z) in the form of a third-order polynomial: This function has been substituted into Equation (16), and the cell potential V cell = V oc − bη 0 withη 0 given by the solution of this equation has been fitted to the perturbed data. Fitting has been performed with the NonlinearFit procedure of Maple . The result is depicted in Figure 4 (the parameters resulted from fitting are a 0 = 0.0162, a 1 = 0.0887, a 2 = 10 −4 and a 3 = 0.691). As can be seen, the quality of fitting is good, and the reconstructed shape ofD(z) is close to that shape used to generate the "experimental" polarization data. A similar procedure has been employed to reconstruct the shape of the exchange current density from the polarization data. Eighteen equidistant points belonging to the "aged" (red) polarization curve in Figure 2b, (which corresponds to the exponential alongz exchange current density) have been generated and perturbed by adding a random 10-mV disturbance to the cell potential. This data was then considered as the "experimental" polarization curve, and the model has been fitted to this data using a third-order polynomial for the parameter, ε: This corresponds to the following fit of the exchange current density: The results of fitting are shown in Figure 5. As can be seen, the quality of polarization curve fitting is high, and thez-shape of i * reconstructed from IV-curve fitting is very close to the exact shape. The coefficients in Equation (23) resulting from fitting are a 1 = 2.55, a 2 = 0.1 and a 3 = 3.51. Note that successful fitting has been obtained by replacing the arcsinh-function in Equation (15) by a logarithm of twice the argument. Overall, these tests show the principal possibility of reconstructing the non-uniform shapes of the CCL parameters along the oxygen channel by accurate fitting of the aged cell polarization curves. Note that the processing of real polarization curves would require that: (i) the curve has been acquired at a large oxygen flow rate and (ii) the curve is IR-corrected. The first requirement is necessary, because calculations with the variable along z oxygen concentration are very time consuming. The second requirement "clears" the polarization curve from the ohmic term, which often nonlinearly depends on the cell current.

Conclusions
We report a quasi-2D model for the polarization curve of a PEM fuel cell with a non-uniform distribution of the CCL parameters along the oxygen channel. The model is applied to simulate the effect of non-uniform CCL aging in PEMFC. The results show that non-uniform distributions of the CCL oxygen diffusivity and exchange current density along the oxygen channel modify the curve in a unique way. The resulting curves differ quite significantly from the curves for the cell with the uniform parameters. We demonstrate numerical reconstruction of the shape of the non-uniform along the channel CCL parameters by fitting the model to the cell polarization data.

Conflicts of Interest
The author declare no conflict of interest.