A quasi three–dimensional model for PEM fuel cell impedance: The effect of local oxygen depletion under the flow field rib

Our previous 1d+1d model for PEM fuel cell impedance is extended to take into account 2d distribution of the oxygen concentration in the gas–diffusion layer under the flow field channel and rib. Fitting of the model impedance to the experimental spectra from the segmented PEMFC shows that in spite of nearly uniform GDL oxygen diffusivity, local oxygen concentration and current density exhibit strong decay under the rib (land). The effect is merely due to larger oxygen transport distance under the rib. The results demonstrate validity of local EIS for measuring 3D distribution of local current in an operating PEM fuel cell.


Introduction
It is well known that water management is an extremely important issue for PEM fuel cell performance [1].In particular, oxygen transport through the porous gas-diffusion layer (GDL) depends on the GDL water content.Large amount of liquid water produced in the ORR could flood the GDL thereby blocking the oxygen transport to adjacent domain of the catalyst layer [2].
In the past two decades, a lot of work has been done to visualize liquid water distribution and dynamics in an operating PEM fuel cell.In operando water imaging is typically done using either X-ray tomography [3][4][5][6][7][8][9][10] or neutron radiography [11][12][13][14][15].However, both techniques require unique and expensive equipment (synchrotron radiation or neutron beam source).Detailed discussion of the techniques for water imaging and key results have been reviewed by Bazylak [16]; more recent works have been discussed by Wang et al [17].
In the context of present work, of particular interest is the distribution of liquid water under the channel and rib (land) GDL areas.Muirhead et al [18] reported X-ray tomography of water distribution in the plane perpendicular to the air channel in the cell equipped with the Sigracet 25 BC GDL.25BC is the carbon paper formed by straight in-plane carbon fibers fixed by PTFE; the paper is coated with MPL. Figure 1 from [18] demonstrates strong non-uniformity of water distribution with high water content under the land (rib) and low water content under the channel.
Similar strong non-uniformity of GDL substrate water content under the rib and channel was reported by Ge et al [19] from X-ray measurements.Ge et al [19] developed also a simple model allowing them to obtain a power in exponential dependence of oxygen transport resistivity of land and channel substrate domains on liquid saturation.Attempt was made to derive the GDL transport resistance by fitting equivalent electric circuit (EEC) to the measured impedance spectra.However, the EEC approach is not reliable, as discussed by Macdonald [20].
Reum et al [21] reported direct measurements of potential distribution along the y-coordinate in figure 1 using the potential probes implanted between the GDL and catalyst layer.Solution of Laplace equation for potential distribution in the GDL with the boundary condition taken from experiment gave the local current distribution.It was shown that under standard operating conditions with relative humidity (RH) between 0.7 and 1.0, local current peaks under the channel.It is worth noting that in [21], the membrane resistance under the channel was found to be 2-3 times lower than under the rib.
Using neutron radiography sequencing technique, Schneider et al [22] studied liquid water transients in the GDL under the channel and rib with simultaneous local current measurements using the segmented anode.They reported strongly nonuniform local current distribution between the channel and rib areas with the peak under the channel.Later Schneider et al [23] developed locally resolved (between rib and channel areas) impedance model.The model included 2d transient oxygen diffusion equation in the GDL and CCL, and zerodimensional transient proton charge conservation equation with the y-coordinate (figure 1) as a parameter.The oxygen diffusion coefficient in the GDL was assumed to be constant, while the local current was dependent on the y-coordinate parametrically, due to parametric variation of the local overpotential along y.The results have shown strong effect of decay of under-rib oxygen concentration and local current on the total cell impedance.Note that the variation of local overpotential along the y-axis would induce proton current along this direction; however, in [23] this current was ignored.
Guan et al [24] developed a pore-network model for oxygen transport through the porous GDL substrate.The substrate structure was obtained from X-ray tomography.The authors demonstrated high sensitivity of the oxygen transport to the presence of water under the land area and relatively weak effect of liquid water under the land (rib) on this transport.Wang et al [17] reported a detailed model of water transport in the cell porous layers.The model employs a high-resolution 3d image of the porous media obtained using X-ray tomography and further enhancement using a deep learning algorithm.Lattice-Boltzmann model was solved in the reconstructed volume and impressive pictures of liquid water transport through the sample with the resolution of less than 1 μm were published.The model shows that in the fabric-based GDL, liquid water is transported to the channel primarily through the network of holes and small cracks, while oxygen is transported through the water-free domains.The model is based on X-ray imaging and it requires extreme computational resources.
Below, we extend our previous 1d+1d impedance model of a PEMFC to take into account 2d transient effects in the plane perpendicular to the channel axis, along the channel/rib coordinate.The oxygen diffusion coefficient in the GDL is allowed to vary along the y-coordinate (figure 1) which seems to be reasonable given the strong non-uniformity of water content.However, fitting of the model to experimental spectra from a segmented cell equipped with the Sigracet 25BC GDL shows that the contribution of oxygen transport under the rib into the cell impedance is marginal.Local oxygen concentration and current density exhibit strong nonuniformity between the channel and rib.The effect is caused by larger oxygen transport length under the rib and it is more pronounced ah higher cell currents.Overall, the results show that analysis of 3D water content in an operating PEM fuel cell with segmented electrodes could be done using relatively inexpensive technique of impedance spectroscopy.Water distribution in the GDL under the channel and land (rib) areas in a PEMFC operating at the limiting current density of about 3 A cm −2 .Red curve-0% RH at the cathode inlet, blue curve-100% RH.Reprinted from Muirhead et al [18] with permission of Elsevier.

Model 2.1. Along the channel and CCL through-plane models
In our previous works, a 1d+1d model of PEMFC impedance has been developed [25][26][27].The model is based on the following idea.One-dimensional oxygen transport equation along the air channel is linked to the onedimensional through-plane mass and charge conservation equations in the GDL and CCL.The link is provided by the volumetric sink term in the channel equation, which models oxygen flux to the GDL due to the oxygen reduction reaction (ORR).
A 1d plug-flow oxygen mass transport equation linked to the oxygen transport in the GDL is depicted in figure 2. This figure also shows the through-plane oxygen and proton conservation equations in the CCL.The channel and CCL sub-models used below do not differ from the previous versions; detailed discussion of equations in figure 2 and boundary conditions can be found in [28].

Oxygen transport in the GDL
Here, we extend the model by taking into account 2d distribution of the oxygen concentration under the channel/rib in the plane perpendicular to the channel axis.As in the previous model, we separate the cell into N segments along the channel coordinate z. Figure 3 illustrates 2d model equations for oxygen transport in the GDL of a single segment: Local oxygen concentration in the channel c h (z) provides a boundary condition at the channel/GDL interface.Under the rib, the normal oxygen flux is zero, and due to problem symmetry it is zero on the left and right sides of the 2d computational domain (figure 3).On the bottom side of the domain, continuity of the normal oxygen flux in the GDL and CCL is applied.The static version of equation ( 1) is solved with the boundary condition at the GDL/CCL interface expressing stoichiometric local normal oxygen flux Here, c b 0 is the static oxygen concentration, l t is the CCL thickness, and j 0 (y) is the local current density, which is assumed to be proportional to the local oxygen concentration at the GDL/CCL interface where k T is the ORR rate constant.Equation (3) is simply the Tafel law under assumption that the ORR overpotential is independent of y.The constant k T is determined from the requirement where L rc is the length of the computational domain along the y-axis and j * (z) is the local current density along the z-coordinate.Quite evidently, where L is the channel length and J is the mean cell current density.Experiments [11,18] indicate that the GDL domain under the rib could be flooded, leading to a nonuniform distribution of the oxygen concentration and local current density j 0 (y) along the y-coordinate (figure 3).To capture this effect, the oxygen diffusion coefficient in the GDL is allowed to vary along the y-axis

Linear equations and impedance
Equation (1) and the equations in figure 2 were transformed using the dimensionless variables Linearization and Fourier-transform of equations in figure 2 is described in [28].The CCL underneath the GDL domain was separated into M sub-segments (yellow rectangles in figure 3) and the CCL problem for the perturbation amplitudes c 1 ˜and 1 h was solved for each sub-segment, taking into account the variation of local static oxygen concentration and current density along the y-coordinate.
Calculations have shown that for all frequencies, the variation of c b 1 ˜along the y-axis is marginal and it can be neglected.Equation (10) thus can be reduced to 1d, through-plane equation which can be easily solved leading to ˜( ), in turn, is used as a boundary condition for the CCL problem on the next iteration.The impedance of the k-th sub-segment was calculated as The impedance of the nth segment Z n ˜(figure 3) was calculated as a sum of parallel sub-segment impedances: Finally, the total segment impedance Z seg,n in the dimension form suitable for fitting was calculated as

Experimental
Local EIS measurements were performed using a segmented cell system developed in Hawaii Natural Energy Institute [29].The setup consists of a custom-built test station, a hardware, current/voltage sensors and PXI data acquisition system (figures 5(a) and (b)).A closed loop Hall sensor from Honeywell (Model CSNN191) was used for the current transducer system which allowed control of up to 10 channels in either a high (standard) or low current mode.The standard current mode enables measurements of segment currents up to 15 A. The current limit of the data acquisition system can be extended to 30 A or more using a counter current technology.The low current mode ensured an accurate segment current measurement up to 400 mA and it was used for cyclic voltammetry and linear sweep voltammetry.Data collection of the voltage and current signals was performed with a National Instrument PXI data acquisition instrument operating under control of HNEI developed LabView programs.The data sampling frequency of the PXI system was 1 MHz, sufficient for measuring simultaneous responses from 10 segments.The segmented fuel cell was operated as a single cell in either voltage or current control mode using HNEIʼs custom-built test station.The total current and power limitations of the system are 240 A and 1.2 kW, respectively.The developed approach mimics real fuel cell operation conditions.The main operating parameters like gas flow rates, humidity, back pressure and overall cell current or voltage are controlled and set by the test station, while measured currents of individual segments are floating.
The segmented cell hardware is built for 100-cm 2 membrane electrode assembly (MEA).The hardware utilizes one set of modified segmented flow field and current collector leaving the other side of the cell nonsegmented.The segmented flow field were divided into 10 segments, where segment 1 is at the inlet, while segment 10 is at the cell outlet.The flow fields on both sides of the cell had the same 10-channel parallel serpentine design and were arranged in co-flow configuration.
Commercially available catalyst coated membrane (CCM) with the active area of 100 cm 2 was used in this work.The catalyst loading for both electrodes was 0.4 mg Pt cm −2 and anode/cathode thickness varied in the range of 10-12 μm.The reinforced membrane had thickness of 16-18 μm.Sigracet 25BC was employed as gas diffusion layer (GDL) for both electrodes.To maintain the optimal compression ratio of the MEA we applied Teflon gaskets of the thickness 125 μm or the anode and cathode.For this work we used segmented GDL and gasket at the cathode electrode as presented at figure 5(c).The electrode area covered by the gasket did not contribute to the produced current which reduced the active MEA area to 76 cm 2 .The in-plane current between segments was estimated to be negligible compared to the current through the segment.
The measurements were performed at 80 °C.The anode and cathode operating conditions were 100 and 50% RH, respectively, 150 kPa back pressure for both electrodes.Hydrogen and air were supplied with stoichiometry of 2 and 4, respectively.The EIS was measured under galvanostatic control of the overall cell with 11 steps/decade.The frequency range was from 10 kHz to 0.05 Hz.The AC perturbation current signal was set in order to get the voltage response amplitude of 10 mV.The cell operating parameters are listed it table 1.

Results and discussion
Equation (2) for the static oxygen concentration was solved on a uniform rectangular computational grid using SSOR method [31].The impedance equation ( 16) was fitted to experimental spectra measured from the ten segments of the segmented cell (section Experimental).Fitting was performed in Python environment using the SciPy least_squares procedure.The AC potential perturbation signal was applied to all segments simultaneously, while the response of segment current was calculated taking into account the transport of oxygen concentration perturbations with the flow along the channel, i.e., the impedance of nth segment was dependent on impedances of all segments located upstream.The details of fitting procedure can be found in [28].
Detailed comparison of experimental and fitted model spectra for the segment #6 at the cell current density of 600 mA cm −2 is shown in figure 6. Figure 7 shows the experimental and fitted Nyquist spectra for all the segments at the cell current of 800 mA cm −2 ; good quality of fitting is seen.).The fitting procedure was rather insensitive to the value of GDL oxygen diffusivity under the rib, meaning that oxygen transport under the rib does not contribute much to the cell performance.This result agrees with the conclusion drawn by Ge et al [19] and Guan et al [24].On the other hand, the value of D b max is nearly uniform along the channel with the minor variation between 0.02 cm 2 s −1 and 0.03 cm 2 s −1 (figure 8(a)).2D distributions of the oxygen concentration c b (x, y) in the segments 1 (air inlet), 5 (middle) and 9 (close to the outlet) together with the distributions of local current density along the y-coordinate are shown in figure 9 for the mean cell current density of 800 mA cm −2 .As can be seen, quite large gradients of the oxygen concentration leads to a strong non-uniformity of the local current density: at the mean current density of about 800 mA cm −2 , local current under the channel is close to 1100 mA cm −2 , while under the rib it drops down to      9).However, at lower current density (200 mA cm −2 ), the oxygen concentration and local current are practically uniform under the channel and rib (figure 10).
As discussed in the Introduction, measurements of liquid water distribution in an operating PEMFC is an expensive and time-consuming task.Figure 9 shows that the effects of local GDL flooding could be studied by means of a much less expensive technique of local electrochemical impedance spectroscopy of segmented cell.

Conclusions
A quasi-3d model for PEM fuel cell impedance is developed.The model takes into account 1d+1d oxygen flow along the channel, 2d oxygen transport through the gas diffusion layer, and 1d oxygen and proton transport through the cathode catalyst layer.Fitting of the model to measured impedance spectra of a segmented PEMFC equipped with the Sigracet 25BC GDL gave the following results.
• The contribution of oxygen transport under the rib to the cell impedance is marginal, which does not allow to reliably determine the oxygen diffusivity of the GDL domain under the rib.Possible local water accumulation under the ribs hence does not affect the effective oxygen transport properties of Sigracet 25BC.
• At the current density of 800 mA cm −2 , oxygen concentration is strongly non-uniform with this concentration under the channel being nearly three times higher than under the rib.Low oxygen concentratiion is due to large oxygen transport distance under the rib.
• At the current density of 800 mA cm −2 , local current density is also strongly non-uniform with almost three times higher value under the channel than under the rib.
• However, at lower cell current density (200 mA cm −2 ) both the oxygen concentration and local current are fairly uniform between the channel and rib.
To summarize, the results show that local EIS is a useful tool for studying 3d local current distribution in an operating PEM fuel cell.

Figure 1 .
Figure 1.Water distribution in the GDL under the channel and land (rib) areas in a PEMFC operating at the limiting current density of about 3 A cm −2 .Red curve-0% RH at the cathode inlet, blue curve-100% RH.Reprinted from Muirhead et al [18] with permission of Elsevier.

Figure 2 .
Figure 2. Equations for the oxygen transport along channel and for the oxygen and proton transport through the CCL.For notations see Nomenclature section.

( 6 )
describes the two-level shape of the oxygen diffusion coefficient with D and the smooth transition zone of the width s between the two values (figure4).Parameters D b min and D b max were determined from fitting (see below).

Figure 3 .Equation ( 1 )
Figure 3. Two-dimensional computational domain (segment) in the plain perpendicular to the channel axis.One-dimensional oxygen transport equation along the channel provides a boundary condition for the 2d oxygen concentration distribution in the xy plane under the rib and channel.Solution of the 2d problem gives the boundary condition for M sub-segments of the CCL; in each sub-segment, the 1d through-plane problem is solved leading to the segment impedance.
taken from solution of the CCL problem in the respective subdomain.The value c 1 b 1

Figure 4 .
Figure 4. Schematic of the oxygen diffusion coefficient shape along the y-axis, equation (6).Parameters D b min and D b max in equation (6) are declared as fitting ones, while s = 0.05 is fixed.
In this work, of largest interest is the effect of 2d oxygen concentration distribution in the GDL.The GDL oxygen diffusivity under the channel D b max and the d-factor in equation D dD b b min max = were declared as fitting parameters.Parameter D b max was determined quite accurately, (figure 8(a)), while the accuracy of d-factor, i.e., of the parameter D b min in all the segments was very poor (figure 8(b)

Figure 5 .
Figure 5. (a) Schematic of the segmented cell system, (b) overall view of the tests station and the segmented cell; (c) segmented flow field.Figure 5(a) is reprinted from [29], figure 5(b) is reprinted from [30], with permission from Elsevier.

Figure 6 .
Figure 6.(a) The experimental (points) and fitted model (open circles) Nyquist spectra of segment #6 for the cell current density of 600 mA cm −2 .The other experimental conditions are collected in table 1.(b) The frequency dependence of the real and imaginary part of the spectrum in (a).

Figure 7 .
Figure 7.The experimental (points) and fitted model (lines) Nyquist spectra of segments 1-10 for the cell current density of 800 mA cm −2 .

Figure 8 .
Figure 8.(a) The distribution of GDL oxygen diffusion coefficient D b max under the channel for the segments 1 to 10 along the channel.(b) d-factor in equation D dD b b min max= for all the segments (see equation (6)).Segment #1 is located at the inlet and segment 10-at the outlet.

Figure 9 .
Figure 9. 2d distribution of the scaled oxygen concentration under the rib and channel in segments 1 (a), 5 (b) and 9 (c) (see figure 3 for the 2d domain).Solid lines: local current density along the rib/channel coordinate y.Mean current density in the cell is 800 mA cm −2 .

Figure 10 .
Figure 10.2d distribution of the scaled oxygen concentration under the rib and channel in segment 5; see figure 3 for the domain.Mean current density in the segment is about 200 mA cm −2 .Solid line: local current density along the rib/channel coordinate y. .

Table 1 .
Pore/cell geometrical and operating parameters.