# Electrochemical Pressure Impedance Spectroscopy for Polymer Electrolyte Membrane Fuel Cells: A Combined Modeling and Experimental Analysis

Lutz Schiffer, Anantrao Vijay Shirsath, Stéphane Raël, Caroline Bonnet, François Lapicque, Wolfgang Bessler

► To cite this version:

Lutz Schiffer, Anantrao Vijay Shirsath, Stéphane Raël, Caroline Bonnet, François Lapicque, et al. Electrochemical Pressure Impedance Spectroscopy for Polymer Electrolyte Membrane Fuel Cells: A Combined Modeling and Experimental Analysis. Journal of The Electrochemical Society, 2022, 169 (3), pp.034503. 10.1149/1945-7111/ac55cd. hal-03614188

HAL Id: hal-03614188

https://hal.univ-lorraine.fr/hal-03614188

Submitted on 9 May 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Electrochemical pressure impedance spectroscopy (EPIS) has recently been developed as a potential diagnosis tool for polymer electrolyte membrane fuel cells (PEMFC). It is based on analyzing the frequency response of the cell voltage with respect to an excitation of the gas-phase pressure. We present here a combined modeling and experimental study of EPIS. A pseudo-two-dimensional PEMFC model was parameterized to a 100 cm<fnref n="2" /> laboratory cell installed in its test bench, and used to reproduce steady-state cell polarization and electrochemical impedance spectra (EIS). Pressure impedance spectra were obtained both in experiment and simulation by applying a harmonic pressure excitation at the cathode outlet. The model shows good agreement with experimental data for current densities ≤ 0.4 A cm⁻<fnref n="2" />. Here it allows a further simulative analysis of observed EPIS features, including the magnitude and shape of spectra. Key findings include a strong influence of the humidifier gas volume on EPIS and a substantial increase in oxygen partial pressure oscillations towards the channel outlet at the resonance frequency. At current densities ≥ 0.8 A cm⁻<fnref n="2" /> the experimental EIS and EPIS data cannot be fully reproduced. This deviation might be associated with the formation and transport of liquid water, which is not included in the model.

© 2022 The Author(s). Published on behalf of The Electrochemical Society by IOP Publishing Limited. This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 License (CC BY, http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse of the work in any medium, provided the original work is properly cited. [DOI: 10.1149/1945-7111/ac55cd]

<figure><img src="image_2.png" /></figure>

Manuscript submitted June 30, 2021; revised manuscript received February 7, 2022. Published March 4, 2022.

List of symbols

A<sub>ch</sub> Cross-sectional area of gas channel, m<fnref n="2" />

a Water activity

a<sub>a,d</sub> Water absorption/desorption rate coefficient, m/s

C<sub>DL</sub><sup>V</sup> Volume-specific double-layer capacitance, F/m³

c<sub>i</sub> Concentration of species i, mol/m³

D<sub>ij</sub><sup>eff</sup> Binary diffusion coefficient of species i and j, m<fnref n="2" />/s

D<sub>ij</sub><sup>ff</sup> Effective binary diffusion coefficient of species i and j, m<fnref n="2" />/s

D<sub>i</sub><sup>K</sup> Knudsen diffusion coefficient of species i, m<fnref n="2" />/s

D<sub>λ</sub> Diffusion coefficient of water in ionomer, m<fnref n="2" />/s

d<sub>particle</sub> Diameter of particle in porous electrode, m

EW Equivalent weight of the dry ionomer, kg/mol

F Faraday's constant, 96485 C/mol

f Volume fraction of water in ionomer

f<sub>1</sub>, f<sub>2</sub> Friction factors of gas channel

ΔH Reaction enthalpy, J/mol

i<sub>cell</sub> Cell current density, A/m<fnref n="2" />

i<sub>elyt</sub> Proton flux, A/m<fnref n="2" />

i<sub>F</sub><sup>V</sup> Volume-specific Faradaic current density, A/m³

i<sub>0</sub><sup>V</sup> Volume-specific exchange current density, A/m³

i<sub>DL</sub><sup>V</sup> Volume-specific current density due to double layer charge/discharge, A/m³

J<sup>conv</sup> Convective molar flux, mol/m<fnref n="2" /> s

J<sub>i</sub><sup>diff</sup> Diffusive molar flux of species i, mol/m<fnref n="2" /> s

J<sub>λ</sub><sup>diff</sup> Diffusive molar flux of water in the ionomer, mol/m<fnref n="2" /> s

J<sub>λ</sub><sup>drag</sup> Molar flux of water dragged in the ionomer, mol/m<fnref n="2" /> s

j<sub>i</sub><sup>diff</sup> Diffusive mass flux of species i, kg/m<fnref n="2" /> s

k<sub>a,d</sub> Water absorption/desorption rate constant, m/s

L<sub>CH</sub> Channel length, m

L<sub>CL</sub> Thickness of the catalyst layers, m

M<sub>i</sub> Molar mass of species i, kg/mol

ṁ Mass flow, kg/s

n<sub>CH</sub> Number of gas channels

P<sub>ch</sub> Channel perimeter, m

p<sub>i</sub> Partial pressure of species i, Pa

P<sub>H₂O</sub><sup>sat</sup> Saturation pressure of water, Pa

p<sub>ref</sub> Reference pressure in Nernst equation, 101325 Pa

R Ideal gas constant, 8.314 J/mol K

R<sub>contact</sub> Specific contact resistance, Ω m<fnref n="2" />

r<sub>pore</sub> Pore radius of porous electrode, m

ΔS Reaction entropy, J/mol K

Ṡ<sub>i</sub> Molar source term of species i, mol/m³ s

Ṡ<sub>λ</sub> Molar source term of water in the ionomer, mol/m³ s

ṡ<sub>i</sub> Mass source term of species i, kg/m³ s

T Temperature, K

T<sub>ref</sub> Reference temperature in kinetic expressions, 353.15 K

t Time, s

V<sub>cell</sub> Cell voltage, V

V<sub>hum</sub> Gas phase volume of humidifier, m³

V<sub>m</sub> Volume per sulfonic acid sides of the ionomer, m³/mol

v Channel flow velocity, m/s

w<sub>CH</sub> Channel width, m

w<sub>rib</sub> Channel rib width, m

x<sub>i</sub> Molar fraction of species i,

Z<sub>V/i</sub> Transfer function: electrochemical impedance, Ω m<fnref n="2" />

Z<sub>V/lp</sub> Transfer function: electrochemical pressure impedance, V/ Pa

Z<sub>p/lp</sub> Transfer function: inlet pressure - outlet pressure,

Z<sub>p(O₂)/lp</sub> Transfer function: average (CCL) oxygen pressure - outlet pressure,

α<sub>A</sub>, α<sub>C</sub> Anodic and cathodic transfer coefficient of half-cell reactions,

ε Porosity of porous electrode,

ε<sub>elyt</sub> Volume fraction of the ionomer,

η Overpotential, V

κ Permeability of the porous electrode, m<fnref n="2" />

λ Water content of the ionomer, -

λ<sup>eq</sup> Water content of the ionomer at equilibrium, -

λ<sub>i</sub> Inlet stoichiometry of species i,

μ Viscosity of gas, kg/m s

ν Frequency, 1/s

ξ Electro-osmotic drag coefficient, -

<footnote n="1">¹Institute of Sustainable Energy Systems (INES), Offenburg University of Applied Sciences, 77652 Offenburg, Germany</footnote>

<footnote n="2">²Reactions and Chemical Engineering Laboratory, CNRS-University of Lorraine, 54000 Nancy, France</footnote>

<footnote>*Electrochemical Society Member.</footnote>

<footnote n="2">²E-mail: lutz.schiffer@hs-offenburg.de; wolfgang.bessler@hs-offenburg.de</footnote>
$\rho_i$ Density of species i, kg/m³

$\rho_{elyt,dry}$ Mass density of the dry ionomer, kg/m³

$\sigma_{elyt}$ Proton conductivity of the ionomer, S/m

$\tau$ Geometric tortuosity of porous electrode, -

$\tau_w$ Wall shear stress, Pa

$\phi$ Electric potential, V

$\phi(Z)$ Phase shift of transfer function Z, °

$\Delta\phi$ Galvani potential, V

$\Delta\phi^{eq}$ Galvani potential at equilibrium, V

$\varphi$ Relative humidity, -

With increasing interest in the hydrogen economy and electric vehicle technologies, fuel cells remain the subject of continuous research and development efforts.¹⁻⁴ Currently, the most prominent type of fuel cell is the polymer electrolyte membrane fuel cell (PEMFC).⁵ The diversity of mass transport processes, coupled with electrochemical reactions, makes the PEMFC a quite complex energy conversion unit. In order to control or optimize these processes, it is crucial to characterize them.

A powerful and widely-used tool to characterize the dynamics of electrochemical systems, such as PEMFCs, is electrochemical impedance spectroscopy (EIS).⁶ By analyzing the electric impedance in response to a small harmonic voltage or current excitation, different processes can be visualized and distinguished from each other in their characteristic frequency range. It is very useful for the characterization of electrochemical processes, for example, charge-transfer kinetics and double layer capacities.⁷ However, the characterization of mass transport processes, commonly assumed to dominate the impedance in the low-frequency impedance region, remains debatable.⁸ Several studies relate the origin of the low-frequency arc to diffusion limitations in the gas diffusion layers or catalyst layers.⁹⁻¹¹ However, Schneider et al.¹²,¹³ experimentally showed that, at least for frequencies below 10 Hz, oxygen concentration oscillations extend into the gas channel. They argue that these oscillations, transported downstream via the gas channel, couple locally with the externally-imposed alternating current (AC). This causes the formation of the low-frequency arc.

In 2010, Niroumand et al. observed periodic fluctuations of the cell voltage with an identical frequency to the fluctuations of the cathode output pressure of the PEMFC at a frequency of 0.14 Hz.¹⁴ In their communication paper, they encouraged researchers to further investigate pressure-induced voltage oscillations at low frequencies in order to develop a diagnostic tool complementary to EIS. This gave rise to the idea of using pressure as a dynamic state variable in addition to current and voltage in frequency response analysis experiments.

For studying pressure dynamics in relation to an electric state variable in the frequency domain, Hartmann et al. introduced in 2014 the term electrochemical pressure impedance spectroscopy (EPIS).¹⁵ In their study of a metal-air battery, they were able to show the sensitivity of EPIS towards transport-related parameters, such as diffusion coefficients. The potential of EPIS to obtain complementary information to standard EIS was shown by Grübl et al. by directly comparing the sensitivities of these two characterization methods to dedicated cell parameters.¹⁶ While EIS was mainly sensitive to electrochemical parameters, EPIS revealed a stronger sensitivity towards transport parameters. The analysis was conducted with a 1D model of a generic half cell comprising a closed gas reservoir.

EPIS was first applied experimentally to a PEMFC in 2016 by Engebretsen et al., who excited the cathode outlet pressure of a 5 cm² cell with a loudspeaker in the frequency range 10 mHz–100 Hz and analyzed the voltage response.¹⁷ Although the spectra were subject to considerable noise, some trends could be observed. At open-circuit voltage (OCV), they reported a ratio of the voltage and pressure amplitudes independent of the excitation frequency and consistent with the Nernst potential variation with pressure. With increasing current density, a substantial increase in the low-frequency magnitude was observed. At constant current density,

the spectra showed a constant magnitude at low frequencies (10 mHz–1 Hz) and a decrease of the magnitude towards zero roughly in the range of 1–100 Hz. The phase shift, on the contrary, was relatively independent of the current density and showed a monotonous decrease with increasing frequency. By varying the inlet humidification they showed that the magnitude increases for dry conditions.

In a more recent study, Shirsath et al. were able to obtain a cleaner signal of the pressure impedance, although only in the frequency range of 1 mHz–1 Hz, by exciting the cathode outlet pressure of a 100 cm² cell with a diaphragm pressure regulator.¹⁸,¹⁹ Qualitatively similar to Engebretsen et al., they observed an increase in the magnitude with current density in the low-frequency region. However, quantitatively these values differ up to one order of magnitude. Furthermore, the change of magnitude with frequency revealed a different qualitative behavior, by showing an increase with frequency up to a maximum slightly below 1 Hz. However, the phase shift was consistent for both setups. Shirsath et al. observed that gas phase volumes present upstream of the cell, particularly inside the gas humidifier, cause a frequency-dependent change of pressure oscillation along the gas channel in terms of amplitude and phase shift.¹⁸,¹⁹

A different approach was followed by Sorrentino et al.²⁰,²¹ They excited the partial pressure of oxygen by periodically adding an oxygen stream at the PEMFC gas inlet and analyzed the relation between the partial pressure and an electric observable such as current or voltage in the frequency domain. This methodology was termed concentration-alternating frequency response analysis (CFRA). They claimed that the spectra are not masked by charge-transfer dynamics and therefore allow a selective analysis of the mass transport processes.

The difference between the different experimental EPIS approaches and their results show the requirement of further studies towards a better understanding of this characterization method. For this purpose, we show here for the first time a combined simulative and experimental analysis of EPIS in a PEMFC. EPIS signals were generated by applying pressure excitation at the cathode outlet. Dynamic models allow an in-depth insight into the dynamic behavior of internal states, in particular concentration and pressure fields. The simulations shown here are based on a dynamic, physics-based, pseudo-two-dimensional model of a single-channel PEMFC. This type of fuel cell model has been broadly used in literature before,²²⁻²⁶ but has not yet been applied to EPIS simulations. We have extended the state-of-the-art by including into the model the gas phase volume in the humidifier installed upstream of the fuel cell in order to account for its previously-observed¹⁸,¹⁹ influence on pressure dynamics. We are thus able to predict EPIS over a wide frequency range (as in Engebretsen et al.¹⁷) while maintaining a high signal quality (as in Shirsath et al.¹⁸,¹⁹). Based on parameter variations and analysis of internal states, the origin of impedance features could be identified.

After presenting the experimental setup and the modeling and simulation methodology, the ability of our model to reproduce experimentally observed EIS and EPIS is shown. There also the origin of specific EPIS features as well as the influence of the gas humidifier are analyzed. Finally, the article is summarized and concluded.

## Methodology

**Experimental setup.**—The experimental setup has been described in detail previously.¹⁸,¹⁹ Briefly, we have used a laboratory-scale 100 cm² single PEMFC (Ulmer Brennstoffzellen-Manufaktur GmbH, UBzM). The gases were guided in counterflow mode through 23 serpentine channels at each flow field. The excitation of pressure at the cathode outlet was realized with a diaphragm-based pressure regulator. All experimental data shown in the present paper were obtained at 55 °C and an average outlet pressure fixed at 116325 Pa in both compartments. The cathode was
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Schematic representation of the 1D+1D modeling domain with indicated mass transport directions through layers cathode/anode gas channel (CCH/ACH), cathode/anode gas diffusion layer (CGDL/AGDL), cathode/anode catalyst layer (CCL/ACL), and polymer electrolyte membrane (PEM).</figcaption></figure>

fed with air at 55% relative humidity and an oxygen stoichiometry of 2.5. The anode was fed with dry hydrogen and a stoichiometry of 1.2.

**Modeling domain.**—The modeling domain of the pseudo-two-dimensional model (P2D or 1D+1D) is shown in Fig. 1. The macroscopic dimension (here: x direction) describes the one-dimensional transport of the gas phase through the gas channels (CH). Although actual flow field geometries (such as a serpentine flow field) are geometrically much more complex, the simple 1D approach still allows the simulation of reactant depletion along the channel. The mesoscopic dimension (here: y direction) describes the one-dimensional transport of mass and charge in the through-plane direction of the membrane electrode assembly (MEA), which is sandwiched between the gas channels. Both dimensions are coupled at their interfaces by setting the boundary fluxes of the MEA as a source term in the gas channel. The MEA itself consists of five different spatially resolved layers: In the middle, separating the cathode and anode, is the polymer electrolyte membrane (PEM). It consists of a proton-conducting ionomer, which also allows the uptake and transfer of water, but is impermeable for gases. Attached to each side of the PEM are the catalyst layers (CL), where the electrochemical reactions occur. In the CL, we assume the coexistence of gas and ionomer phases as well as an electron-conducting phase. The outer layers of the MEA are the gas diffusion layers (GDL), in which the gas phase is transported through the pores of an electron-conducting matrix. The geometric parameters (thicknesses, porosity, tortuosity, particle and pore diameters, electrode composition) of CL, GDL and channels are assumed identical for anode and cathode.

It was observed before¹⁸,¹⁹ that the humidifier used in the experimental PEMFC setup for humidifying the inflow gas streams has a strong influence on the EPIS signal. Therefore, a 0D model of a humidifier gas volume upstream of the cathode channel was included.

The PEMFC model captures major transport processes in the along-the-channel direction and through-the-MEA dimension, both of which are relevant for EPIS. However, the model neglects the presence of liquid water inside the cell; water in the pores of the electrodes and the gas channels is assumed to be always gaseous. The model is therefore strictly valid only under non-condensing conditions (sufficiently low inflow humidification and/or sufficiently low current densities). The investigation of the impact of liquid-water presence and transport on EPIS signal will be the subject of future investigations. The model is isothermal.

In the following Sections, the governing equations of electrochemistry and transport are described, followed by the parametrization of the model.

**Governing equations.**—Electrochemistry.—The half-cell reactions taking place at anode and cathode are modeled as single-step

reactions. The reaction in the anode catalyst layer (ACL) is the hydrogen oxidation reaction (HOR),

$$H_2(g) \rightleftharpoons 2 H^+(elyt) + 2 e^-(elde). \quad [1]$$

Here, “elyt” and “elde” stand for electrolyte (the ionomer) and electrode (electron conductor), respectively. The reaction in the cathode catalyst layer (CCL) is the oxygen reduction reaction (ORR), where we assume water to form dissolved in the ionomer phase,²⁷,²⁸

$$\frac{1}{2}O_2(g) + 2 H^+(elyt) + 2 e^-(elde) \rightleftharpoons H_2O(elyt). \quad [2]$$

The local kinetics of both half-cell reactions are described by relating the volume-specific Faradaic current density $i_F^V$ to the overpotential $\eta$ with the Butler-Volmer equation,²⁷

$$i_F^V = i_0^V \left( \exp\left[\frac{\alpha_A F}{RT}\eta\right] - \exp\left[-\frac{\alpha_C F}{RT}\eta\right] \right), \quad [3]$$

where, $i_0^V$ is the volumetric exchange current density. Instead of using a symmetry factor which is strictly valid only for single-step reactions, we use the anodic and cathodic transfer coefficients $\alpha_A$ and $\alpha_C$. These and the other symbols are defined in the symbol list. The local volumetric Faradaic current is defined as positive for an oxidation and negative for a reduction reaction. The overpotential driving the electrochemical reactions is the difference between the local Galvani potential $\Delta\phi$ and the Galvani potential at equilibrium $\Delta\phi^{eq}$,²⁹

$$\eta = \Delta\phi - \Delta\phi^{\text{eq}}. \quad [4]$$

The Galvani potential is the potential difference between electron-conducting and proton-conducting phase,²⁹

$$\Delta\phi = \phi_{\text{elde}} - \phi_{\text{elyt}}. \quad [5]$$

The equilibrium potential difference in the respective catalyst layer is given by the Nernst equation,²⁷

$$\Delta\phi^{\text{eq}} = \begin{cases} -\frac{\Delta H_{\text{ORR}} - T \Delta S_{\text{ORR}}}{2F} + \frac{RT}{2F} \ln\left(\sqrt{\frac{p_{O_2}}{p_{\text{ref}}}}\right), & \text{CCL} \\ \\ \frac{\Delta H_{\text{HOR}} - T \Delta S_{\text{HOR}}}{2F} - \frac{RT}{2F} \ln\left(\frac{p_{H_2}}{p_{\text{ref}}}\right), & \text{ACL.} \end{cases} \quad [6]$$

The volume-specific current density due to charge and discharge of the double layer caused by potential transients is calculated as,²⁹

$$i_{\text{DL}}^{V} = C_{\text{DL}}^{V} \frac{\partial(\Delta\phi)}{\partial t}, \quad [7]$$

where, $C_{\text{DL}}^V$ is the volume-specific double-layer capacitance.

**Charge transport.**—The charge is transported as electrons in the electron-conducting phase and as protons in the ionomer phase. CL and GDL are assumed to ideally conduct the electrons, resulting in a spatially constant potential of the electron-conducting phase.²⁹ However, the limited conductivity of the ionomer $\sigma_{\text{elyt}}$ is a significant cause for voltage loss. The potential gradient of the ionomer phase is related to the proton flux by Ohm’s law,²⁹

$$\nabla\phi_{\text{elyt}} = -\frac{i_{\text{elyt}}}{\sigma_{\text{elyt}}}. \quad [8]$$

The charge continuity equation relates the local area-specific proton flux $i_{\text{elyt}}$ to the local current production rates due to Faradaic current and electric double layer charge and discharge,²⁹
$$ \nabla i_{\text{elyt}} = i_{\text{F}}^{\text{V}} + i_{\text{dl}}^{\text{V}}. \quad [9] $$

The macroscopic observables are the cell voltage $V_{\text{cell}}$ and the cell current density $i_{\text{cell}}$. Accounting for contact resistances $R_{\text{contact}}$, the cell voltage is given as

$$ V_{\text{cell}} = \phi_{\text{elde,C}} - \phi_{\text{elde,A}} - i_{\text{cell}} \cdot R_{\text{contact}}. \quad [10] $$

The cathode potential $\phi_{\text{elde,C}}$ is set to zero as reference point.²⁹ The cell current follows from the integration of the local current densities over the complete cathode,

$$ i_{\text{cell}} = -\frac{1}{L_{\text{CH}}} \int_{x=0}^{L_{\text{CH}}} \int_{y=0}^{L_{\text{CL}}} (i_{\text{F}}^{\text{V}} + i_{\text{DL}}^{\text{V}}) dy dx. \quad [11] $$

Boundary conditions for the transport Eq. 8 are given in Table I.

*Gas-phase transport in the channels.*—The mass transport occurs in the gas channels and pores of the MEA as gas, assumed as an ideal gas, and in the ionomer as dissolved water. Within the gas channels we assume that the transport can be described sufficiently as one-dimensional transport along the channel,³⁰ the predominant transport mechanism is convection. The bulk continuity equation is

$$ \frac{\partial \rho}{\partial t} = -\nabla(\rho v) + \sum_{i=1}^{N} \dot{s}_i. \quad [12] $$

Here and in the following, $i$ is the index for species. The exchange of gas-phase species between the channel and the adjacent GDL pores is considered with the source term $\dot{s}_i$, which is connected to the gas flux of the porous layer boundary,

$$ \dot{s}_i = -\frac{w_{\text{CH}} + w_{\text{rib}}}{A_{\text{CH}}} J_i^{\text{CH}|\text{GDL}} M_i, \quad [13] $$

where the boundary flux is defined to be positive in direction of the GDL. Note that, because of the 1D description of the channel, the mass transport over the channel/electrode boundary appears as source term of the channel continuity equation, not as boundary condition.

We additionally consider free molecular diffusion, so that the continuity equation for species $i$ is²⁹,³¹

$$ \frac{\partial \rho_i}{\partial t} = -\nabla(\rho_i v) - \nabla j_i^{\text{diff}} + \dot{s}_i. \quad [14] $$

By assuming a 1D flow, neglecting the gravity forces, and assuming that the viscous forces are dominated by wall shear stress, the momentum conservation equation simplifies to³¹

$$ \frac{\partial(\rho v)}{\partial t} = -\nabla(\rho v^2) - \nabla p - \frac{P_{\text{ch}}}{A_{\text{ch}}} \tau_w. \quad [15] $$

The diffusive flux is calculated with an analytical solution of the Maxwell-Stefan equations for multicomponent diffusion by using binary diffusion coefficients³¹ according to

$$ \nabla x_i = - \sum_{j=1}^{N} \frac{1}{c D_{ij}} \left( \frac{x_j j_i^{\text{diff}}}{M_i} - \frac{x_i j_j^{\text{diff}}}{M_j} \right). \quad [16] $$

In order to close the set of equations the sum of diffusive mass fluxes is set to zero,

$$ \sum_{i} j_{i}^{\text{diff}} = 0, \quad [17] $$

which implies that the velocity $v$ in the gas channels describes the mass average velocity.³¹ Boundary conditions are given in Table I.

*Gas-phase transport in the porous layers.*—The transport mechanisms of gas in the pores of the GDL and CL are convection and diffusion, therefore, the species continuity is given as

$$ \frac{\partial(\epsilon c_i)}{\partial t} = -\nabla(J^{\text{conv}}x_i) - \nabla J_i^{\text{diff}} + \dot{S}_i. \quad [18] $$

Note that, different to the channel continuity equation, we use the molar concentration (instead of the mass density) as conservative variable.²⁹,³⁰ The source term (needed in the CL, but not in the GDL) represents the production or consumption of species due to reactions and, in case of water, the phase change by absorption/desorption into respectively out of the ionomer.²⁷,²⁹

<table><caption>**Table I. Boundary conditions for the governing equations of charge and mass transport.**</caption><thead><tr><th>Electrolyte charge transport</th><th>i<sub>elyt</sub><sup>CGDL|CCL</sup> = i<sub>elyt</sub><sup>ACL|AGDL</sup> = 0</th></tr></thead><tbody><tr><td>Gas transport channel</td><td>v<sub>CCH</sub><sup>in</sup> and v<sub>ACH</sub><sup>in</sup> see Eq. 35<br/>p<sub>CCH</sub><sup>out</sup> see Table II, varies harmonically during EPIS<br/>p<sub>ACH</sub><sup>out</sup> see Table II</td></tr><tr><td>Gas transport electrode</td><td>x<sub>O2,CCH</sub><sup>in</sup> = x<sub>O2,CCH</sub><sup>in,dry</sup> (1 - x<sub>H2O,CCH</sub><sup>in</sup>)<br/>x<sub>N2,CCH</sub><sup>in</sup> = x<sub>N2,CCH</sub><sup>in,dry</sup> (1 - x<sub>H2O,CCH</sub><sup>in</sup>)<br/>x<sub>H2O,CCH</sub><sup>in</sup> = q<sub>CCH</sub><sup>in</sup> p<sub>H2O</sub><sup>sat</sup> / p<sub>CCH</sub><sup>in</sup><br/>x<sub>H2,ACH</sub><sup>in</sup> = x<sub>H2,ACH</sub><sup>in,dry</sup> (1 - x<sub>H2O,ACH</sub><sup>in</sup>)<br/>x<sub>H2O,ACH</sub><sup>in</sup> = q<sub>ACH</sub><sup>in</sup> p<sub>H2O</sub><sup>sat</sup> / p<sub>ACH</sub><sup>in</sup></td></tr><tr><td>Dissolved water transport ionomer</td><td>J<sub>i</sub><sup>diff,CCL|PEM</sup> = J<sub>i</sub><sup>diff,PEM|ACL</sup> = 0<br/>J<sub>i</sub><sup>conv,CCL|PEM</sup> = J<sub>i</sub><sup>conv,PEM|ACL</sup> = 0<br/>c<sub>i</sub><sup>CCH|CGDL</sup> = c<sub>CCH,i</sub><br/>c<sub>i</sub><sup>ACH|AGDL</sup> = c<sub>ACH,i</sub><br/>p<sup>CGDL|CCH</sup> = p<sub>CCH</sub><br/>p<sup>AGDL|ACH</sup> = p<sub>ACH</sub><br/>J<sub>λ</sub><sup>CGDL|CCL</sup> = J<sub>λ</sub><sup>ACL|AGDL</sup> = 0</td></tr></tbody></table>
The convective transport through the porous electrodes, driven by the pressure gradient, is given by Darcy's law,³¹

$$v = -\frac{\kappa}{\mu} \nabla p. \qquad [19]$$

We convert the flow velocity to molar fluxes to be used in the species continuity equation,

$$J^{\text{conv}} = c_v. \qquad [20]$$

Note that Darcy's law is commonly used to model pressure-driven convective flow in fuel cell electrodes, and was applied before either in a mass density setting,³² in a molar concentration setting²⁹,³⁰,³³ (as here), or in a mixed setting.³⁴ There are also models that neglect Darcy flow and use diffusion as only transport mechanism.²⁷ It is beyond the scope of the present article to further investigate this difference.

Similar to the gas transport in the gas channels, the diffusive flux through the pores is calculated by the use of the Maxwell-Stefan equations,³¹

$$\nabla x_i = - \sum_{j=1}^{N} \frac{1}{c D_{ij}^{\text{eff}}} (x_j J_i^{\text{diff}} - x_i J_j^{\text{diff}}). \qquad [21]$$

A difference to the diffusion in the gas channels is the reference of the diffusive flux to the average molar velocity instead of the average mass velocity. The closing equation of the Maxwell-Stefan equations therefore is,³¹

$$\sum_j J_i^{\text{diff}} = 0. \qquad [22]$$

The effective diffusion coefficients $D_{ij}^{\text{eff}}$ are calculated based on Bosanquet's formula to account for free molecular and Knudsen diffusion. By calculating averaged diffusion coefficients, the symmetry of the effective binary diffusion coefficients is ensured,²⁹

$$D_{ij}^{\text{eff}} = D_{ji}^{\text{eff}} = \frac{\epsilon}{\tau^2} \frac{1}{2} \left( \frac{1}{1/D_i^K + 1/D_{ij}} + \frac{1}{1/D_j^K + 1/D_{ji}} \right). \qquad [23]$$

Here, $\tau$ is the geometric tortuosity and $\tau^2$ the tortuosity factor.³⁵ Knudsen diffusion describes transport when molecules are colliding more often with the wall than with other molecules. The diffusion coefficient, therefore, is a function of the pore size $r_{\text{pore}}$,³¹

$$D_i^K = \frac{2}{3} r_{\text{pore}} \left( \frac{8 RT}{\pi M_i} \right)^{\frac{1}{2}}. \qquad [24]$$

In the following, the source terms for the species continuity equation of the catalyst layers are defined. In the ACL, the hydrogen consumption due to the HOR is given as

$$\dot{S}_{\text{H}_2} = -\frac{i_F^V}{2F}. \qquad [25]$$

Similarly, in the CCL, the oxygen consumption due to the ORR is given as

$$\dot{S}_{\text{O}_2} = \frac{i_F^V}{4F}. \qquad [26]$$

The rate of water exchanging between gaseous and dissolved form in the ionomer via absorption/desorption is depending on the saturation of the ionomer and is modeled as³⁶,³⁷

$$\dot{S}_{\text{H}_2\text{O}} = -\frac{k_{\text{a,d}}}{L_{\text{CL}} V_{\text{m}}} (\lambda_{\text{eq}} - \lambda), \qquad [27]$$

where, $V_m = EW/\rho_{\text{elyt,dry}}$ is the volume per sulfonic acid sides of the ionomer. Absorption occurs if the water content of the ionomer is smaller than the equilibrium water content $\lambda^{\text{eq}}$. Boundary conditions are given in Table I.

*Water transport in the ionomer.*—Inside the ionomer, the flux of dissolved water is modeled based on the work of Springer et al.³⁸ The two transport mechanisms are diffusion and electro-osmotic drag. We formulate the water continuity equation using the water uptake $\lambda$ as a state variable, defined as the number of water molecules per acid group, according to

$$\frac{\partial \lambda}{\partial t} \frac{\epsilon_{\text{elyt}}}{V_{\text{m}}} = -\nabla J_{\lambda}^{\text{diff}} - \nabla J_{\lambda}^{\text{drag}} + \dot{S}_{\lambda}, \qquad [28]$$

where, $\epsilon_{\text{elyt}}$ is the volume fraction of the ionomer. The dissolved water diffuses through the ionomer due to a gradient of water content,

$$J_{\lambda}^{\text{diff}} = -\frac{D_{\lambda}}{V_{\text{m}}} \nabla \lambda. \qquad [29]$$

The protons drag water molecules from anode to cathode, which is expressed as

$$J_{\lambda}^{\text{drag}} = \frac{\xi}{F} i_{\text{elyt}}. \qquad [30]$$

In the CL, the water can enter or leave the ionomer via absorption or desorption, and water is generated in the cathode catalyst layer due to the ORR, which is accounted for by the source term $\dot{S}_{\lambda}$.²⁷ We model it as³⁶

$$\dot{S}_{\lambda} = \begin{cases} \frac{k_{\text{a,d}}}{L_{\text{CL}} V_{\text{m}}} (\lambda^{\text{eq}} - \lambda) - \frac{i_F^V}{2F}, & \text{CCL} \\ \frac{k_{\text{a,d}}}{L_{\text{CL}} V_{\text{m}}} (\lambda^{\text{eq}} - \lambda), & \text{ACL}. \end{cases} \qquad [31]$$

Boundary conditions are given in Table I.

*Humidifier.*—PEMFCs typically include gas humidifiers upstream of the cell, which are needed to keep the membrane sufficiently hydrated. We have observed¹⁸ that the gas volume of the humidifier strongly influences the pressure dynamics along the gas channel. Therefore, we include it explicitly in the model. The mass conservation in the humidifier gas volume, represented as continuously stirred tank reactor, is given as

$$\frac{\partial \rho_{\text{hum}}}{\partial t} = \frac{n_{\text{CH}}(\dot{m}_{\text{feed}} - \dot{m}_{\text{CH}}^{\text{in}})}{V_{\text{hum}}}, \qquad [32]$$

where, the feed $\dot{m}_{\text{feed}}$ is the mass influx of the humidified gas. The number of flow-field channels $n_{\text{CH}}$ has to be considered because the fuel cell model describes only a single channel, while the humidifier supplies all channels. The mass flow rate out of the humidifier is assumed to equal the mass flow rate $\dot{m}_{\text{CH}}^{\text{in}}$ into the fuel cell channels. It can be expressed as

$$\dot{m}_{\text{CH}}^{\text{in}} = v_{\text{CH}} A_{\text{CH}} \rho_{\text{CH}}^{\text{in}}. \qquad [33]$$

By assuming constant temperature and no pressure losses in gas lines between humidifier and fuel cell, the gas density in the humidifier $\rho_{\text{hum}}$ equals the density of the gas at the channel inlet $\rho_{\text{CH}}^{\text{in}}$,
$$ \rho_{\text{hum}} = \rho_{\text{CH}}^{\text{in}}. \quad [34] $$

Combining Eqs. 32–34 results in the following expression for the cathode inlet velocity,

$$ v_{\text{CH}}^{\text{in}} = \frac{\left( \dot{m}_{\text{feed}} - \frac{V_{\text{hum}}}{n_{\text{CH}}} \frac{\partial \rho_{\text{CH}}^{\text{in}}}{\partial t} \right)}{A_{\text{ch}} \rho_{\text{CH}}^{\text{in}}}, \quad [35] $$

which we use as a boundary condition for the channel flow model, Eq. 15.

The mass flow of the humidified gas $\dot{m}_{\text{feed}}$ is regulated according to the respective stoichiometry $\lambda_i$ ($i$ = O₂, H₂),

$$ \dot{m}_{\text{feed}} = \frac{\lambda_i}{x_{i,\text{CH}}^{\text{in}}} \frac{i A_{\text{FC}}}{z_i F n_{\text{CH}}} \bar{M}_{\text{CH}}^{\text{in}} \quad [36] $$

with

$$ A_{\text{FC}} = (w_{\text{CH}} + w_{\text{rib}}) L_{\text{CH}} n_{\text{CH}}. \quad [37] $$

As the present system only considers humidification of the cathode gas stream, the equation above simplifies for the anode with $V_{\text{hum}} = 0$. We assumed the anode gas feed to have a remaining humidity of $\phi_{\text{ACH}}^{\text{in}} = 0.05$.

**Parameterization.**—Electrochemistry.—In order to solve the Butler-Volmer equation, the thermodynamic and kinetic coefficients need to be parameterized for both half-cell reactions. The thermodynamic properties needed to calculate the equilibrium potential for both reactions are given in Table II. For the volumetric exchange current density of the ORR, the following expression is used, in which the pressure and temperature dependence is based on Ref. 39, while the volumetric exchange current density at reference conditions has been fitted to our experiments,

$$ i_0^V = 1.55 \cdot 10^6 \frac{A}{m^3} \left( \frac{p_{O_2}}{p_{\text{ref}}} \right)^{0.54} \exp \left( \frac{67 \frac{\text{kJ}}{\text{mol}}}{R} \left( \frac{1}{T_{\text{ref}}} - \frac{1}{T} \right) \right) \quad [38] $$

This expression shows that a change in oxygen concentration due to mass transport limitations is affecting the reaction kinetics and consequently, the cell voltage. The cathodic transfer coefficient of the ORR has been fitted to the experiments as $\alpha_C = 1.434$. The anodic transfer coefficient was calculated according to the original formulation of the BV-equation in Ref. 27 as $\alpha_A = 2 \cdot (1 - \frac{\alpha_C}{2}) = 0.566$.

The volumetric exchange current density of the HOR is calculated according to the measurements made by Durst et al.²⁷,⁴² with an anode platinum loading used in Ref. 27,

$$ i_0^V = 2.7 \cdot 10^3 \frac{A}{m_{\text{Pt}}^2} \cdot 10^7 \frac{m_{\text{Pt}}^2}{m^3} \exp \left( \frac{16 \frac{\text{kJ}}{\text{mol}}}{R} \left( \frac{1}{T_{\text{ref}}} - \frac{1}{T} \right) \right). \quad [39] $$

The transfer coefficients $\alpha_A$ and $\alpha_C$ of the HOR are set to 1.²⁷ The volume-specific double-layer capacitance $C_{\text{DL}}^V$ at the cathode was fitted to the subsequently shown EIS results and assumed to be equal at the anode. The value is given in Table II.

**Charge transport.**—As we are assuming ideal electronic conductivity in the electrodes, only the proton conductivity of the

ionomer needs to be parametrized. The proton conductivity is calculated according to Springer et al.,³⁸ who fitted the conductivity as a function of the water content to experiments conducted at 30 °C. The activation energy was derived from an additional measurement of a fully hydrated membrane at 80 °C and is assumed to be valid for all water contents. In combination with the Bruggeman correction to account for the volume fraction of ionomer in the catalyst layers,⁴³ the conductivity is given as

$$ \sigma_{\text{elyt}} = \epsilon_{\text{elyt}}^{1.5} (0.5139\lambda - 0.326) \frac{S}{m} \exp \left( 1268 \left( \frac{1}{303} - \frac{1}{T} \right) \right). \quad [40] $$

For $\lambda < 1$, the conductivity is assumed constant at the value of $\lambda = 1$.

**Mass transport.**—For the gas channels, the wall shear stress of the momentum conservation equation is expressed as a function of the velocity according to

$$ \tau_w = \frac{1}{8} \frac{P_{\text{CH}}}{A_{\text{CH}}} \mu v f_1 + \frac{1}{2} \rho v^2 f_2, \quad [41] $$

where, the first term describes the pressure loss for laminar flow in a square duct, and the second term the pressure losses due to additional flow hindrance such as turns of the serpentine flow field. Both factors $f_1$ and $f_2$ were fitted to experimental pressure losses measured at 20 °C and 50 °C and are given in Table II.

In order to calculate the convective gas transport in the electrodes, the permeability of the porous media is needed. It is calculated with the Kozeny-Carman equation, which expresses the permeability as a function of geometric tortuosity, porosity, and particle diameter according to³⁰

$$ \kappa = \frac{\epsilon^3 d_{\text{particle}}^2}{72\tau(1 - \epsilon)^2}. \quad [42] $$

The binary gas diffusion coefficients $D_{ij}$ are taken from Fuller et al.⁴⁴

The diffusivity of dissolved water in the ionomer is calculated based on an expression derived by Motupally et al.⁴⁵ They fitted the diffusivity to the measured values of Zawodzinski et al. at 30 °C⁴⁶ and added a temperature correction obtained by Yeo and Eisenberg.⁴⁷ With consideration of the volume fraction of the ionomer in the catalyst layers by the Bruggeman correction,⁴³ the diffusivity is

$$ D_{\lambda} = \begin{cases} 3.1 \cdot 10^{-7} \frac{m^2}{s} \cdot \epsilon_{\text{elyt}}^{1.5} \lambda (-1 + \exp(0.28\lambda)) \exp\left(-\frac{2436 \text{ K}}{T}\right), & \text{for } \lambda \le 3 \\ 4.17 \cdot 10^{-8} \frac{m^2}{s} \cdot \epsilon_{\text{elyt}}^{1.5} \lambda (1 + 161 \cdot \exp(-\lambda)) \exp\left(-\frac{2436 \text{ K}}{T}\right), & \text{for } \lambda > 3. \end{cases} \quad [43] $$

The electro-osmotic drag coefficient, defining the number of water molecules dragged per proton, is assumed to be linearly proportional to the water content, as suggested by Springer et al.,⁵⁸

$$ \xi = \frac{2.5\lambda}{22}. \quad [44] $$

The water transfer coefficients for absorption and desorption are assumed to increase proportionally with the volume fraction of water,³⁶

$$ k_{\text{a,d}} = a_{\text{a,d}} f \exp \left[ \frac{20 \frac{\text{kJ}}{\text{mol}}}{R} \left( \frac{1}{T_{\text{ref}}} - \frac{1}{T} \right) \right], \quad [45] $$
<table><caption>**Table II.** Model parameters, taken from literature (indicated by reference), fitted to experimental data (indicated by asterisk *), and obtained from experimental setup.</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td><strong>Geometry MEA</strong></td><td></td></tr><tr><td>Thickness L {GDL, CL, PEM}</td><td>{212.5, 10, 50}·10<sup>-6</sup> m</td></tr><tr><td>Porosity ϵ {GDL, CL}</td><td>{0.77, 0.4}<sup>27</sup></td></tr><tr><td>Tortuosity τ {GDL, CL}</td><td>{1.26, 1.26}</td></tr><tr><td>Volume fraction ionomer ϵ<sub>elyt</sub> {CL, PEM}</td><td>{0.3, 1.0}<sup>27</sup></td></tr><tr><td>Particle diameter d<sub>particle</sub> {GDL, CL}</td><td>{5.4, 5.4}·10<sup>-7</sup> m</td></tr><tr><td>Pore radius r<sub>pore</sub> {GDL, CL}</td><td>{2.7, 2.7}·10<sup>-7</sup> m</td></tr><tr><td><strong>Geometry gas channel</strong></td><td></td></tr><tr><td>Channel length L<sub>CH</sub></td><td>0.304 m</td></tr><tr><td>Channel cross-section area A<sub>CH</sub></td><td>3.08·10<sup>-7</sup> m<sup>2</sup></td></tr><tr><td>Channel width w<sub>CH</sub></td><td>7·10<sup>-4</sup> m</td></tr><tr><td>Channel rib width w<sub>rib</sub></td><td>7.3·10<sup>-4</sup> m</td></tr><tr><td>Channel wetted perimeter P<sub>CH</sub></td><td>2.28·10<sup>-3</sup> m</td></tr><tr><td>Number of channels n<sub>CH</sub></td><td>23</td></tr><tr><td>Gas phase volume in cathode humidifier V<sub>hum</sub></td><td>850 ml</td></tr><tr><td><strong>Ionomer</strong></td><td></td></tr><tr><td>Equivalent weight EW</td><td>1.02 kg/equiv<sup>27</sup></td></tr><tr><td>Density of dry ionomer ρ<sub>elyt,dry</sub></td><td>1.97·10<sup>-3</sup> kg/m<sup>327</sup></td></tr><tr><td><strong>Electrochemistry and transport</strong></td><td></td></tr><tr><td>Reaction enthalpy ΔH<sub>ORR</sub>, ΔH<sub>HOR</sub></td><td>-285.83 kJ/mol<sup>40</sup>, 0kJ/mol</td></tr><tr><td>Reaction entropy ΔS<sub>ORR</sub>, ΔH<sub>HOR</sub></td><td>-163.3 J/mol K<sup>41</sup>, 0.104 J/mol K<sup>41</sup></td></tr><tr><td>Exchange current density i<sub>0</sub><sup>V</sup> (p<sub>ref</sub>, T<sub>ref</sub>) of HOR</td><td>2.7·10<sup>-3</sup> A/m<sub>pt</sub><sup>-2</sup>·10<sup>7</sup> m<sub>pt</sub><sup>2</sup>/m<sup>327,42</sup></td></tr><tr><td>Exchange current density i<sub>0</sub><sup>V</sup> (p<sub>ref</sub>, T<sub>ref</sub>) of ORR</td><td>1.55·10<sup>6</sup> A/m<sup>3*</sup></td></tr><tr><td>Transfer coefficients α<sub>A</sub>, α<sub>C</sub> of HOR</td><td>1<sup>27</sup>, 1<sup>27</sup></td></tr><tr><td>Transfer coefficients α<sub>A</sub>, α<sub>C</sub> of ORR</td><td>1.434, 0.566*</td></tr><tr><td>Double layer capacity C<sub>DL</sub><sup>V</sup></td><td>4.1·10<sup>7</sup> F/m<sup>3*</sup></td></tr><tr><td>Contact resistance R<sub>contact</sub></td><td>6.3·10<sup>-6</sup> Ω m<sup>2*</sup></td></tr><tr><td>Friction factor f<sub>1</sub></td><td>18.5*</td></tr><tr><td>Friction factor f<sub>2</sub></td><td>0.021*</td></tr><tr><td><strong>Operation conditions</strong></td><td></td></tr><tr><td>Temperature T</td><td>328.15 K</td></tr><tr><td>Outlet pressure p<sub>CH</sub><sup>out</sup></td><td>116325 Pa</td></tr><tr><td>Inlet relative humidity ρ<sub>ACH</sub><sup>in</sup>, ρ<sub>CCH</sub><sup>in</sup></td><td>0.05, 0.55</td></tr><tr><td>Inlet mixture (dry) anode x<sub>H2</sub></td><td>1</td></tr><tr><td>Inlet mixture (dry) cathode x<sub>O2</sub>, x<sub>N2</sub></td><td>0.21, 0.79</td></tr><tr><td>Inlet stoichiometry λ<sub>H2</sub>, λ<sub>O2</sub></td><td>1.2, 2.5</td></tr><tr><td><strong>Constants</strong></td><td></td></tr><tr><td>Reference temperature T<sub>ref</sub></td><td>353.15 K</td></tr><tr><td>Reference pressure p<sub>ref</sub></td><td>101325 Pa</td></tr><tr><td>Ideal gas constant R</td><td>8.314 J/mol K</td></tr><tr><td>Faraday constant F</td><td>96485 C/mol</td></tr></tbody></table>

with the kinetic absorption parameter $a_a = 3.53 \cdot 10^{-5} \frac{\text{m}}{\text{s}}$ and the
desorption parameter $a_d = 1.42 \cdot 10^{-4} \frac{\text{m}}{\text{s}}$. The volume fraction of
water in the ionomer is a function of the water content and calculated
according to

$$f = \frac{\lambda V_w}{\lambda V_w + V_m}, \qquad [46]$$

with the volume per sulfonic acid sides of the ionomer $V_m$ and the
molar volume of water $V_w = M_{H_2O}/\rho_{H_2O}$.

For the equilibrium water content of the ionomer, a polynomial is
taken from Springer et al.,³⁸ who fitted the polynomial to experi-
ments at a temperature of 30 °C,

$$\lambda^{\text{eq}} = \begin{cases} 0.043 + 17.81a - 39.85a^2 + 36.0a^3, & \text{for } a \le 1, \\ 14 + 1.4(a - 1), & \text{for } 1 < a \le 3 \end{cases}, \quad [47]$$

with the water activity $a = p_{H_2O}/p_{H_2O}^{\text{sat}}(T)$. Similar to Springer et al., who used the expression at 80 °C, we assume it not to change up to our operating temperature of 55 °C. If the partial pressure of water exceeds the saturation pressure, the sorption is extended by a linear term allowing a maximum water content of 16.8.

**Impedance spectroscopy simulations.**—In the present study, impedance spectroscopy simulations were conducted with different types of excitation and response variables. Simulations of “standard” EIS were performed galvanostatically using a step excitation of the fuel cell current density with a 10% increase in the current load within 0.1 µs, followed by a 10⁵ s relaxation period of the potential.⁴⁸ The electrochemical impedance was calculated as the relation between the Fourier-transformed voltage response and current excitation,

$$Z_{V/\dot{V}}(\nu) = \frac{\mathcal{F}\{\Delta V_{\text{cell}}(t)\}}{\mathcal{F}\{\Delta i_{\text{cell}}(t)\}}. \quad [48]$$
In order to reproduce the experiments, the inlet mass flow rate was kept constant during the step excitation, according to the stoichio-metry of the steady-state current density. Also, the outlet pressure was kept constant.

Pressure impedance spectroscopy simulations were performed under galvanostatic control while the outlet pressure of the cathode gas channel was excited harmonically with an amplitude of 100 Pa. For any given frequency, the harmonic excitation was repeated until the change of electrochemical pressure impedance (defined subse- quently) magnitude was less than 1% compared to the previous excitation wave. Simulations were carried out for excitation fre- quencies between 1 mHz and 100 Hz with four samples per decade. In order to test if the applied outlet pressure excitation evokes a linear voltage response, the pressure amplitude was reduced from 100 Pa to 50 Pa. The change of amplitude and phase of the electrochemical pressure impedance in the whole spectrum was less than 0.8%, so that the excitation amplitude of 100 Pa can be assumed to be sufficient small.

In the following, different transfer functions used to describe the cell response to the pressure excitation are defined. The so-called electrochemical pressure impedance was calculated as the relation between the Fourier-transformed voltage response and pressure excitation,

$$Z_{V/p}(\nu) = \frac{\mathcal{F}\{\Delta V_{\text{cell}}(t)\}}{\mathcal{F}\{\Delta P_{\text{cch}}^{\text{out}}(t)\}}. \quad [49]$$

Further transfer functions of internal state responses to the pressure excitation were calculated as defined in the following. The relation of the Fourier-transformed cathode inlet pressure response to the pressure excitation was obtained as

$$Z_{p/p}(\nu) = -\frac{\mathcal{F}\{\Delta P_{\text{cch}}^{\text{in}}(t)\}}{\mathcal{F}\{\Delta P_{\text{cch}}^{\text{out}}(t)\}}. \quad [50]$$

The relation of the Fourier-transformed average oxygen pressure response in the cathode catalyst layer to the pressure excitation was obtained as

$$Z_{p(O_2)/p}(\nu) = -\frac{\mathcal{F}\{\Delta \bar{P}_{\text{ccl,O}_2}(t)\}}{\mathcal{F}\{\Delta P_{\text{cch}}^{\text{out}}(t)\}}. \quad [51]$$

In order to allow a clear distinction of these different impedances in the text and the figures, we use the symbols $Z_{V/i}$, $Z_{V/p}$, $Z_{p/p}$ and $Z_{p(O_2)/p}$ throughout the remainder of this article.

**Numerical implementation.**—The fuel cell model was imple- mented in the in-house multiphysics simulation package DENIS (“detailed electrochemistry and numerical impedance simulation”).<sup>23,29</sup> Mathematically, the model is represented by a system of partial differential equations and algebraic equations. The spatial derivatives were discretized using a finite-volume scheme with five control volumes in the x scale and 29 control volumes in the y scale. The resulting differential-algebraic equation system was time-integrated using the adaptive, semi-implicit time-step solver LIMEX (version 2.4).<sup>49,50</sup>

DENIS is a C/C<sup>++</sup> code with interfaces to MATLAB. All simulations, including the impedance simulations and the required post-processing and visualization, were carried out in MATLAB (version R2019a).

## Results and Discussion

In this section, we first present the ability of the model to reproduce experimentally observed steady-state cell polarization, gas channel pressure losses, EIS, and pressure impedance spectra. For low current loads, we further show further EPIS simulations in

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Experimental (solid circles) and simulated (solid line) polarization curve at 55 °C, 101325 Pa outlet pressure, with humidified air at the cathode (55% RH, λ<sub>O₂</sub> = 2.5) and dry hydrogen at the anode (5% RH, λ<sub>H₂</sub> = 1.2, cf. Table II).</figcaption></figure>

which we study partial pressure oscillations inside the cell and the influence of the humidifier gas volume.

**Polarization curve.**—Experimental and simulated polarization curves of the analyzed fuel cell are shown in Fig. 2, noting that the experiments were only carried out up to a current density of 1 A cm<sup>-2</sup> due to safety settings of the bench. The simulation shows a good qualitative agreement with the experiment. A difference becomes visible above a current density of 0.4 A cm<sup>-2</sup>, where the slope of the simulated curve begins to decrease continuously, while the experimental slope remains nearly constant. The experiments show a significantly lower OCV than the simulation, which predicts the Nernst equilibrium voltage. There are several possible reasons for the origin of the low OCV, which have been discussed by Vilekar and Datta.<sup>51</sup> These authors conclude that the gas transfer between the electrodes is the most likely reason, which is not included in the present model. Overall, the results confirm the ability of the model to reproduce steady-state current-voltage relations.

**Electrochemical impedance spectroscopy.**—In order to validate the dynamic current-voltage behavior and to identify unknown model parameters, in particular, the double-layer capacitance $C_{DL}^V$ and the contact resistance $R_{contact}$, we have simulated EIS for different current densities. Results are shown in comparison to experimental data in Fig. 3 in Nyquist and Bode representations. The experimental spectra reveal two overlapping arcs in the Nyquist plot, which are more pronounced at lower current densities. With increasing current density, the high-frequency arc (HFA) decreases and is progressively masked by the increasing low-frequency arc (LFA). It is well-accepted in literature<sup>52</sup> that the HFA can be related to the dynamics of the ORR and the LFA to mass transfer processes. For the model, we have fitted the double-layer capacitance $C_{DL}^V$ of the CCL such that the characteristic time of the HFA matches the experiments for a current density of 0.2 A cm<sup>-2</sup>. The impedance spectrum at this current shows a qualitatively good agreement with the experiments, although the LFA is less pronounced. With increasing current density, the experimental and simulated spectra diverge. While the LFA in the experiments increases with increasing current density and masks the HFA, the LFA of the simulation (order of 1 Hz) gets reduced, and another arc at even lower frequencies (order of 10 mHz) evolves, indicating an additional dynamic process. The contact resistance $R_{contact}$ of the model was fitted
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Experimental (solid circles) and simulated (solid lines) electrochemical impedance spectra $Z_{V/i}$ in (a) Nyquist and (b), (c) Bode representation for different current densities between 0.2 and 1.0 A cm⁻². The experimental DC resistance (gradient of the polarization curve) is plotted on the real axis of the Nyquist plot (crosses). The operating conditions are the same as for the polarization experiments shown in Fig. 2.</figcaption></figure>

such that the simulated ohmic resistance, visible in the Nyquist plot as high-frequency intercept with the abscissa, shows the same value as the experiments for a current density of 0.2 A cm⁻². In the experiments, the ohmic resistance is independent of current density. In contrast, the simulated ohmic resistance increases with current density. As the fitted contact resistance is constant, the change in the ohmic resistance has to be caused by an increase in the membrane resistance, resulting from a decreasing membrane hydration. The inductive behavior of the experimental data above 1 kHz is expected to be caused by the connecting wires.⁶

The total AC resistance of the cell, shown as the low-frequency limit of the real part of the impedance, is significantly higher in the experiments compared to the simulations (up to a factor of almost two). For further analysis, we have added the slope of the experimental steady-state polarization curve, also referred to as direct current (DC) resistance, as crosses to Fig. 3a). These values should be equal to the low-frequency impedance measurements, which is clearly not the case for the experiments. The difference between AC and DC resistance has been previously observed and related to the change of stoichiometry (constant mass flow according to the stoichiometry of the average current density) during EIS measurements.⁵³,⁵⁴ However, as we carried both EIS and polarization simulations with the same mass flow conditions as in the experiments, this effect cannot explain the difference between simulations and experiments.

To conclude, with the present model, the dynamic current-voltage relationship can be reproduced qualitatively well for a current density of 0.2 A cm⁻². With increasing current density, the qualitative agreement decreases, probably caused by the different dynamics of mass transport-related processes.

**Gas channel pressure losses.**—The accurate description of pressure losses in the flow field is a requirement for reliable interpretation of EPIS. Figure 4 shows experimental and simulated pressure drops along the cathode channels as a function of air flow rate for different temperatures and outlet pressures at OCV. The two factors $f_1, f_2$ used to describe the wall shear stress (Eq. 41) were fitted such that the model can reproduce the experiments accurately over the complete range of investigated conditions.

**Pressure impedance spectroscopy.**—In the following, we aim to show the ability of the model to reproduce experimentally observed cell response to cathode outlet pressure excitation. As shown above, we have defined two representations of transfer functions, that is, the ratio between cathode inlet and outlet pressure oscillation $Z_{p/p}$, and

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Experimental (solid circles) and simulated (solid lines) pressure drop in the cathode channels as function of air flow rate (total flow rate of all $n_{CH}$ = 23 channels) for different temperatures and outlet pressures at OCV.</figcaption></figure>

the ratio between cell voltage and cathode outlet pressure oscillation $Z_{V/p}$ (also referred to as EPIS).

$Z_{p/p}$ results are shown for different current densities in the four upper panels of Fig. 5. They allow studying the propagation of pressure excitation waves along the cathode gas channel. Experimental data are available in the frequency range of 1 mHz–1 Hz, while we extended the simulations up to 100 Hz. The simulations show a very good agreement with the experiments over the complete frequency range. The amplitude ratio $|Z_{p/p}|$ (panel a, e) approaches a low-frequency limit of almost unity and decreases with increasing frequency down to zero. At low frequencies, this means that the gas pressure oscillates nearly homogeneously throughout the whole channel length. The low-frequency limit is below unity for the following reason. The channel is subject to a constant mass inflow rate; therefore, for increasing outlet pressure, the flow velocity decreases; and as pressure loss depends on velocity (cf. Eq. 41), the inlet pressure shows a smaller increase than the outlet pressure. The effect increases with increasing mass flow rate and therefore, as the mass flow rate is lambda-controlled, with increasing current density. This effect becomes visible in the simulations (panel a) but not in the experiments (panel e).
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Bode plot of simulated (solid lines) and experimental (solid lines with circles) pressure impedance spectra at different current densities. (a), (b), (c), (f) Transfer function $Z_{p/p}$, ratio between cathode inlet and outlet pressure oscillation. (c), (d), (g), (h) Transfer function $Z_{V/p}$, electrochemical pressure impedance. The data at 0.4 A/cm² have been published before.¹⁸</figcaption></figure>

For frequencies larger than ca. 10 mHz, the results show that the pressure oscillation at the inlet gets increasingly dampened until the pressure excitation signal cannot reach the gas channel inlet anymore above 10 Hz or so. The phase shift (panel b, f) exhibits a low-frequency limit of 0°, decreases to -90° around 1 Hz, and slightly increases again towards higher frequencies. It can be interpreted as the delay of the pressure oscillation between the inlet and outlet. At low frequencies, the inlet pressure follows the outlet excitation without delay (low phase shift), but with increasing frequency, this pressure increasingly lags behind the outlet excitation. It should be noted that the behavior of $Z_{p/p}$ is strongly related to the gas humidifier upstream of the inlet, as discussed below. The results shown in the four upper panels of Fig. 5 demonstrate that the frequency-dependent propagation of the pressure excitation wave along the gas channel can be reproduced accurately with the simulation by incorporating the humidifier gas phase volume in the model.

These results form the basis for further analyzing the effect of pressure excitation on the cell voltage, that is, EPIS. A comparison of experimental and simulated EPIS at different current densities is shown in the four lower panels of Fig. 5. First the variations of the absolute value $|Z_{V/p}|$ (panel c, g) are discussed. The low-frequency limit, which can also be referred to as quasi-static electrochemical pressure impedance, shows small values in the order of 1 µV/Pa. The quasi-static electrochemical pressure impedance increases with increasing current density, which is more pronounced in the experiments than in the simulations. Engebretsen et al.¹⁷ reported the same trend but observed values up to one order of magnitude

higher (Engebretsen: 40 µV/Pa at 1.2 A cm⁻², this work 1.9 µV/Pa at 1.0 A cm⁻²). However, they reported a steady-state change of voltage with outlet pressure of 2.4 µV/Pa at 1.6 A/cm², which is in the same order of magnitude as our observation. This means that the quasi-static value might not have been reached in the EPIS measurements of Engebretsen et al.

With increasing frequency, our experimentally-observed $|Z_{V/p}|$ remains constant up to ca. 10 mHz, followed by an increase up to a maximum for frequencies slightly below 1 Hz. The origin of this maximum will be further discussed below. With the present model, the experiments can be reproduced qualitatively at low current density. For current densities $\ge$ 0.8 A cm⁻², the simulated and experimental data diverge. In particular, opposite to the experiments, the simulations show a decrease in $|Z_{V/p}|$ with increasing frequency. This is consistent with the findings of EIS results (cf. Fig. 3), indicating again that the model does not represent all processes occurring at high current density. At further increasing frequency, the simulated $|Z_{V/p}|$ shows a decrease towards zero. This behavior is also indicated by the present experimental data and matches with the observed decrease of pressure impedance in the region of 1–100 Hz by Engebretsen et al.¹⁷

Figures 5d and 5h show the variations of phase shift $\phi(Z_{V/p})$ with the excitation frequency. In the low-frequency limit (quasi-static behavior), the simulations show a phase shift of 0°. This means that a quasi-static increase in the cathode backpressure leads to a proportional increase in the cell voltage. Zhang et al. pointed out that this is both due to increasing equilibrium potential and
decreasing overpotentials.⁵⁵ For increasing frequencies, the experi-
mental data show a strictly monotonic decrease from 0° to about
-200° at the highest measured frequency of 1 Hz, with a slightly
stronger decrease for lower currents. The simulations are in good
quantitative agreement with the experiments. The data show that the
maximum of |Z<sub>V/p</sub>| occurs at the same frequency where
φ(Z<sub>V/p</sub>) = -180°. At further increasing frequencies, the simulated
phase shift continues to decrease to below -400° at the highest
simulated frequency of 100 Hz. This is similar to what was observed
by Engebretsen et al.¹⁷ Physically, this means that the dynamic
response of the voltage is increasingly delayed with respect to
pressure excitation for increasing frequencies.

**Partial pressure analysis.**—In order to further interpret the observed EPIS behavior, we have analyzed the partial pressure oscillations in the CCL. These data are not accessible experimentally but are readily available from the internal states of the model.

Figure 6 shows partial pressures of O₂, N₂, and H₂O as well as total pressure inside the CCL (averaged over the CCL thickness) at three different positions along the channel (inlet, center, and outlet) as a function of time during one excitation period for a current density of 0.2 A cm⁻². The panels on the left side were taken at 1 mHz (quasi-static excitation frequency) and the panels on the right side at 0.56 Hz (resonance frequency, maximum |Z<sub>V/p</sub>|). In the quasi-static case (Figs. 6a–6c), the total pressure curves are similar for all three positions along the channel, meaning that pressure oscillation is independent of the *x* direction and corresponds to the excitation amplitude (100 Pa). The change in partial pressure corresponds to the product of total pressure change and molar fraction, except that the partial pressure of water is constant at the inlet (upper panel) due to the boundary condition of constant relative humidity. The behavior at 0.56 Hz (Figs. 6d–6f) is considerably more complex. At the channel inlet (upper panel), the water partial pressure is almost constant, and the pressure of the other species show only low amplitudes. Towards the channel center (middle panel), the amplitudes strongly increase (note the different *y* scale in the graphs with 0.56 Hz in comparison with those for 1 mHz). At the channel outlet (lower panel), the total pressure oscillation is still given by the external pressure excitation (100 Pa). However, the oxygen pressure oscillation is strongly amplified towards the channel outlet, the amplitude being more than three times higher than the excitation

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Bode plot of simulated pressure impedance spectra for a transfer function Z<sub>p(O<sub>2</sub>)/p</sub> (ratio between CCL-averaged oxygen partial pressure and outlet pressure) at different current densities.</figcaption></figure>

amplitude and the phase being reversed to around -180°. In the
following we refer to this observation as oxygen pressure resonance.
The nearly total phase inversion is not observed for the other species
(nitrogen, water vapor) whose fluctuations follow the excitation
signal. As the total pressure oscillation is a fixed boundary condition
at the outlet, the oscillation of the other species must compensate the
oxygen oscillation. A similar amplification of oxygen pressure
oscillation along the air channel has been reported previously in
the context of EIS measurements experimentally¹² and
theoretically.⁵⁴

As the cell voltage is related to oxygen partial pressure via the
Nernst equation (Eq. 6) and via the exchange current density of the
ORR (Eq. 38), the observed maximum in EPIS is likely related to the
oxygen pressure resonance. In order to further investigate this
relationship, we have calculated the additional transfer function
Z<sub>p(O<sub>2</sub>)/p</sub> according to Eq. 51. It describes the ratio between the

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Species partial pressures, averaged over the CCL thickness, during one excitation period *T* at gas channel inlet, center, and outlet at 0.2 A cm⁻². The amplitude of the outlet pressure excitation was 100 Pa.</figcaption></figure>
oxygen pressure (averaged over the complete CCL, both through the CCL thickness and along the channel length) and the outlet pressure. The resulting spectra are shown in Fig. 7. They show a single feature with a resonance frequency of 0.56 Hz at 0.2 A cm⁻², which decreases at increasing current density. By comparing the spectra to EPIS shown above (Fig. 5), it becomes clear that the characteristic behavior of the EPIS is mainly following the change in oxygen pressure in both amplitude and phase. One characteristic difference is visible in the low-frequency region of the magnitude: while |*Z*<sub>*V*<sub>/p</sub></sub>| is increasing with current density, the |*Z*<sub>*p*(*O*₂)</sub>| is nearly independent of the cell load. Towards low frequencies, |*Z*<sub>*p*(*O*₂)</sub>| is close to 0.15, corresponding to the average molar fraction of oxygen in the CCL at the considered current density (0.2 A cm⁻²). Therefore, it can be concluded that the EPIS resonance is caused by the resonance in oxygen partial pressure.

**Humidifier gas volume.**—This Section is dedicated to the influence of the gas phase volume of the cathode humidifier on pressure impedance spectra. Figure 8 shows spectra of the transfer functions Z<sub>p/p</sub> and Z<sub>V/p</sub> for humidifier gas volumes between 0 ml (no gas volume in the humidifier) and 1000 ml. Note that the simulations in the previous sections were carried out with a volume of 850 ml, corresponding to the experimental setup.

The transfer function Z<sub>p/p</sub> describing the ratio between the inlet and outlet pressure oscillation at 0.2 A cm<sup>-2</sup> (Figs. 8a, 8b) exhibits a strong dependence on the gas phase volume. With increasing volume, both the dampening |Z<sub>p/p</sub>| and delay φ(Z<sub>p/p</sub>) of the pressure excitation shift towards lower frequencies. Without considering the humidifier (V<sub>hum</sub> = 0), the excitation amplitude remains constant along the channel throughout nearly the complete investigated frequency range. At the same time, the profiles of the phase shift vs frequency changes: for high volumes, the curves show a minimum near -90°, while for low volumes, they continue decreasing with increasing frequency.

The electrochemical pressure impedance Z<sub>V/p</sub> at 0.2 A cm<sup>-2</sup> (Figs. 8c, 8d) is also strongly influenced by the humidifier gas volume. Smaller gas volumes lead to lower maxima of the magnitude |Z<sub>V/p</sub>| and a shift towards higher frequencies. For volumes below 10 ml, even the qualitative behavior of the magnitude changes such that the maximum completely vanishes, meaning that there is no more oxygen pressure resonance (Fig. 8c). The phase shift φ(Z<sub>V/p</sub>) also shows a qualitative change. Volumes above 2 ml reveal a strict monotonic decrease with increasing frequency (Fig. 8d), leading to the same phase shift of smaller than -400° at 100 Hz, whereas volumes of 2 ml and less show a relatively smaller non-

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Bode plot of simulated pressure impedance spectra for the transfer function Z<sub>V/p</sub> (electrochemical pressure impedance) at low current densities without the humidifier. The electrochemical double layer (DL) is switched either on (solid lines) or off (dashed lines).</figcaption></figure>

monotonic decrease of phase shift up to a value of less than -100° at 100 Hz.

The simulations show that the cell response to the outlet pressure excitation is strongly connected to the humidifier gas volume upstream of the cell. For the transfer function Z<sub>p/p</sub>, this has previously been pointed out.<sup>18,19</sup> The EPIS resonance only occurs if a certain gas phase volume in the humidifier is present. It is therefore interesting to investigate the case without humidifier gas phase volume in more detail. Results are shown in Fig. 9 for 0.2 and 0.4 A cm<sup>-2</sup>. The peak of |Z<sub>V/p</sub>| at around 0.5 Hz is not present anymore, instead |Z<sub>V/p</sub>| decreases monotonously over the complete investigated frequency range. The phase angle φ(Z<sub>V/p</sub>) remains close to 0° for frequencies up to around 1 Hz before it drops. A difference between 0.2 and 0.4 A cm<sup>-2</sup> is mostly visible at frequencies below around 10 mHz.

In order to investigate the influence of the electrochemical double layer on EPIS, additional simulations were carried out in which the double-layer capacitance was set to zero. This was done in order to investigate the influence of a purely electrochemical process (the

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Bode plot of simulated pressure impedance spectra for different cathode humidifier gas volumes at 0.2 A cm⁻². (a), (b) Transfer function Z<sub>p/p</sub>, ratio between cathode inlet and outlet pressure oscillation. (c), (d) Transfer function Z<sub>V/p</sub>, electrochemical pressure impedance.</figcaption></figure>
double layer charging and discharging) on the pressure impedance. In this case, the pressure impedance reveals two inflection points in the magnitude and a local minimum in the phase shift, indicating two time-dependent processes. The influence of the double layer becomes visible above 1 Hz, leading to a further decrease of magnitude and continuous decrease of phase shift and therefore masking of the higher frequency process. This is noteworthy as EPIS has been described before to depend mainly on transport processes,¹⁶ while the present results show an evident influence of the double layer, too.

## Conclusions

In the present study, a combined simulative and experimental analysis of EPIS has been presented. For current densities $\le 0.4$ A cm⁻², the model shows good agreement with experimentally observed steady-state cell polarization and EIS data. By incorporating a 0D model for the gas volume in the cathode humidifier, experimentally observed cell response to pressure excitation, especially electrochemical pressure impedance, could be successfully reproduced by the model. The experimentally observed dampening and delaying of the pressure excitation wave along the gas channel with increasing frequency were shown to originate from the humidifier gas volume, which turned out to be crucial for the interpretation of EPIS. Both experiments and simulations show a maximum of electrochemical pressure impedance magnitude with respect to the excitation frequency. The maximum is associated with a substantial increase in oxygen pressure oscillation towards the channel outlet, locally up to three times the excitation pressure amplitude.

The model was not able to fully reproduce EIS and EPIS response for current densities $\ge 0.8$ A cm⁻². Here the experimental EIS show an additional low-frequency feature, and the EPIS shows an inverse behavior (magnitude increase) to the model (magnitude decrease). A possible origin of these effects is the formation and transport of liquid water, which is not considered in the present model. If this should be indeed the case, then EPIS is quite sensitive to liquid-water formation. The direction of future studies should be both the model extension to formation and flow of liquid water in the porous structures and a more detailed interpretation of EPIS behavior, for example, through model reduction techniques.

## Acknowledgments

The authors acknowledge funding from the German Research Foundation (Deutsche Forschungsgemeinschaft, DFG) (grant BE 3819/6-1) and from ANR (National Agency for Research in France) (project 17-CE05-0031), both in the framework of the EPISTEL project.

## ORCID

Lutz Schiffer https://orcid.org/0000-0002-0187-5031

Anantrao Vijay Shirsath https://orcid.org/0000-0003-0827-5677

Stéphane Raël https://orcid.org/0000-0001-6642-1607

Caroline Bonnet https://orcid.org/0000-0001-6931-5121

François Lapicque https://orcid.org/0000-0001-9834-7341

Wolfgang G. Bessler https://orcid.org/0000-0001-8037-9046

## References

1. D. Banham, J. Zou, S. Mukerjee, Z. Liu, D. Yang, Y. Zhang, Y. Peng, and A. Dong, *J. Power Sources*, **490**, 229515 (2021).

2. Y. Song, C. Zhang, C.-Y. Ling, M. Han, R.-Y. Yong, D. Sun, and J. Chen, *Int. J. Hydrogen Energy*, **45**, 29832 (2020).

3. Y. Wang, K. S. Chen, J. Mishler, S. C. Cho, and X. C. Adroher, *Appl. Energy*, **88**, 981 (2011).

4. H. Yuan, H. Dai, X. Wei, and P. Ming, *J. Power Sources*, **468**, 228376 (2020).

5. E4tech, *The Fuel Cell Industry Review*, (2020), (https://fuelcellindustryreview.com/).

6. S. M. Rezaei Niya and M. Hoorfar, *J. Power Sources*, **240**, 281 (2013).

7. A. Lasia, *Electrochemical Impedance Spectroscopy and its Applications* (Springer, New York, NY) (2014).

8. S. J. Andreasen, J. R. Vang, and S. K. Kær, *Int. J. Hydrogen Energy*, **36**, 9815 (2011).

9. T. E. Springer, T. A. Zawodzinski, M. S. Wilson, and S. Gottesfeld, *J. Electrochem. Soc.*, **143**, 587 (1996).

10. J. Zhang, Y. Tang, C. Song, and J. Zhang, *J. Power Sources*, **172**, 163 (2007).

11. N. H. Jalani, M. Ramani, K. Ohlsson, S. Buelte, G. Pacifico, R. Pollard, R. Staudt, and R. Datta, *J. Power Sources*, **160**, 1096 (2006).

12. I. A. Schneider, S. A. Freunberger, D. Kramer, A. Wokaun, and G. G. Scherer, *J. Electrochem. Soc.*, **154**, B383 (2007).

13. I. A. Schneider, D. Kramer, A. Wokaun, and G. G. Scherer, *Electrochem. Commun.*, **154**, B770 (2007).

14. A. M. Niroumand, W. Mérida, M. Eikerling, and M. Saif, *Electrochem. Commun.*, **12**, 122 (2010).

15. P. Hartmann, D. Grübl, H. Sommer, J. Janek, W. G. Bessler, and P. Adelhelm, *J. Phys. Chem. C*, **118**, 1461 (2014).

16. D. Grübl, J. Janek, and W. G. Bessler, *J. Electrochem. Soc.*, **163**, A599 (2016).

17. E. Engebretsen, T. J. Mason, P. R. Shearing, G. Hinds, and D. J. L. Brett, *Electrochem. Commun.*, **75**, 60 (2017).

18. A. V. Shirsath, S. Raël, C. Bonnet, L. Schiffer, W. G. Bessler, and F. Lapicque, *Curr. Opin. Electrochem.*, **20**, 82 (2020).

19. A. V. Shirsath, S. Raël, C. Bonnet, and F. Lapicque, *Electrochim. Acta*, **363**, 137157 (2020).

20. A. Sorrentino, T. Vidakovic-Koch, and K. Sundmacher, *J. Power Sources*, **412**, 331 (2019).

21. A. Sorrentino, T. Vidakovic-Koch, R. Hanke-Rauschenbach, and K. Sundmacher, *Electrochim. Acta*, **243**, 53 (2017).

22. M. Eschenbach, R. Coulon, A. A. Franco, J. Kallo, and W. G. Bessler, *Solid State Ionics*, **192**, 615 (2011).

23. J. P. Neidhardt, D. N. Fronczek, T. Jahnke, T. Danner, B. Horstmann, and W. G. Bessler, *J. Electrochem. Soc.*, **159**, A1528 (2012).

24. V. Gurau, H. Liu, and S. Kakaç, *AICHE J.*, **44**, 2410 (1998).

25. A. Z. Weber and J. Newman, *J. Electrochem. Soc.*, **151**, A326 (2004).

26. C. Bao and W. G. Bessler, *J. Power Sources*, **275**, 922 (2015).

27. R. Vetter and J. O. Schumacher, *Comput. Phys. Commun.*, **234**, 223 (2019).

28. H. Wu, X. Li, and P. Berg, *Electrochim. Acta*, **54**, 6913 (2009).

29. W. G. Bessler, S. Gewies, and M. Vogler, *Electrochim. Acta*, **53**, 1782 (2007).

30. H. Zhu, R. J. Kee, V. M. Janardhanan, O. Deutschmann, and D. G. Goodwin, *J. Electrochem. Soc.*, **152**, A2427 (2005).

31. R. B. Bird, W. E. Stewart, and E. N. Lightfoot, *Transport phenomena* (Wiley, New York, NY) (2007) 2nd ed., p. 80.

32. L. M. Pant, S. K. Mitra, and M. Secanell, *Int. J. Heat Mass Transfer*, **58**, 70 (2013).

33. J. H. Nam and D. H. Jeon, *Electrochim. Acta*, **51**, 3446 (2006).

34. H. Zhu and R. J. Kee, *J. Electrochem. Soc.*, **155**, B715 (2008).

35. N. Epstein, *Chem. Eng. Sci.*, **44**, 777 (1989).

36. S. Ge, X. Li, B. Yi, and I.-M. Hsing, *J. Electrochem. Soc.*, **152**, A1149 (2005).

37. H. Wu, *Mathematical Modeling of Transient Transport Phenomena in PEM Fuel Cells* (UWSpace, Ontario, Canada) (2009), http://hdl.handle.net/10012/4862.

38. T. E. Springer, T. A. Zawodzinski, and S. Gottesfeld, *J. Electrochem. Soc.*, **138**, 2334 (1991).

39. K. C. Neyerlin, W. Gu, J. Jorne, and H. A. Gasteiger, *J. Electrochem. Soc.*, **153**, A1955 (2006).

40. M. W. Chase, *JANAF thermochemical tables* (American Inst. of Physics, New York, NY) (1986).

41. M. J. Lampinen and M. Fomino, *J. Electrochem. Soc.*, **140**, 3537 (1993).

42. J. Durst, C. Simon, F. Hasché, and H. A. Gasteiger, *J. Electrochem. Soc.*, **162**, F190 (2015).

43. J. Dujc, A. Forner-Cuenca, P. Marmet, M. Cochet, R. Vetter, J. O. Schumacher, and P. Boillat, *J. Electrochem. En. Conv. Stor.*, **15**, 021001 (2018).

44. E. N. Fuller, P. D. Schettler, and J. C. Giddings, *Ind. Eng. Chem.*, **58**, 18 (1966).

45. S. Motupally, A. J. Becker, and J. W. Weidner, *J. Electrochem. Soc.*, **147**, 3171 (2000).

46. T. A. Zawodzinski, M. Neeman, L. O. Sillerud, and S. Gottesfeld, *J. Phys. Chem.*, **95**, 6040 (1991).

47. S. C. Yeo and A. Eisenberg, *J. Appl. Polym. Sci.*, **21**, 875 (1977).

48. W. G. Bessler, *J. Electrochem. Soc.*, **154**, B1186 (2007).

49. R. Ehrig, U. Nowak, L. Oeverdieck, and P. Deuflhard, "Advanced Extrapolation Methods for Large Scale Differential Algebraic Problems." *High performance scientific and engineering computing*, ed. H.-J. Bungartz (Springer, Berlin) Lecture notes in Computational Science and Engineering, p. 233 (1999).

50. P. Deuflhard, E. Hairer, and J. Zugck, *Numer. Math.*, **51**, 501 (1987).

51. S. A. Vilekar and R. Datta, *J. Power Sources*, **195**, 2241 (2010).

52. X. Yuan, H. Wang, J. Colin Sun, and J. Zhang, *Int. J. Hydrogen Energy*, **32**, 4365 (2007).

53. M. Chandesris, C. Robin, M. Gerard, and Y. Bultel, *Electrochim. Acta*, **180**, 581 (2015).

54. G. Maranzana, J. Mainka, O. Lottin, J. Dillet, A. Lamibrac, A. Thomas, and S. Didierjean, *Electrochim. Acta*, **83**, 13 (2012).

55. J. Zhang, C. Song, J. Zhang, R. Baker, and L. Zhang, *J. Electroanal. Chem.*, **688**, 130 (2013).