## Abstract

This work presents a pseudo-two-dimensional proton exchange membrane fuel cell (PEMFC) model incorporating Nafion ionomer structure–property relationships in the cathode catalyst layer (CL) to capture and explain losses at low Pt loading. Structural data from neutron reflectometry and thin film Nafion conductivity measurements predict variations in the oxygen diffusion coefficient and ionic conductivity with changing CL ionomer thickness and Pt loading. By including these structure–property relationships, predicted polarization curves agree closely with previously published experimental data from cells with Pt loadings between 0.025 and 0.2 mg/cm^{2}. Results demonstrate that structure–property relationships based on physically measurable ionomer and CL properties provide a feasible interpretation of PEMFC CL phenomena for a range of Pt loadings and help explain previously unaccounted-for losses at low Pt. Results also show that simulations must account for surface species coverage variations in order to properly capture the kinetic losses. Finally, results suggest that an increase in ionomer thickness surrounding the C/Pt surfaces may lead to improved cell performance due to improved ionic conductivity.

## 1 Introduction

Proton exchange membrane fuel cells (PEMFCs) provide clean and efficient conversion of chemical energy into electrical energy. Fuel cell electric vehicles (FCEVs) are proposed as a competitor to battery electric vehicles (B-EVs) for decarbonized transportation. FCEVs have an advantage in categories such as range, weight, and refueling times when compared to B-EVs [1]. Although FCEVs are currently available, challenges in infrastructure, performance, durability, and cost hinder wide-spread adoption [2–6].

The PEMFC performance, cost, and durability are all closely tied to the Pt catalyst used in the fuel cell electrodes. Efforts to lower the amount of Pt loading [7–9] or to replace the catalyst altogether [10,11] have shown unacceptable cell performance or low durability. The U.S. Department of Energy (DOE) has set a goal to lower the PEMFC Pt loading to 0.125 mg/cm^{2} or 0.15 mg/kW. However, current state-of-the-art PEMFCs with loadings of roughly 0.13 mg/cm^{2} [12] suffer from performance losses. Electrodes with lower Pt loadings exhibit large transport losses that are not yet sufficiently explained in the literature [13]. Some researchers propose resistances at the ionomer/gas interface and ionomer/Pt interface as a source of these resistances [14,15], although other work suggests limited impact for the former [16]. Although incorporating these resistances provides reasonable agreement to data for existing models, their cause is not yet fully understood.

Multiple PEMFC limiting phenomena occur in the cathode catalyst layer (CL), where nano-thin sheets of ionomer (<10 nm thick) coat carbon-supported Pt nanoparticles. Voltage losses in the CL are attributed to water transport, flooding, electro-osmotic drag, resistances associated with the Pt surface, and transport through the nano-thin ionomer [17–19]. Because of the difficulty in experimentally measuring transport properties in films at this scale, modelers commonly use parameters taken from bulk ionomer materials. However, recent literature demonstrates significant differences between nano-thin and bulk-scale ionomers such as Nafion. Limited direct measurements of the ionic conductivity have been reported [20,21], and the thin film oxygen diffusion coefficients are largely unreported due to experimental limitations. That said, recent studies provide significant insight into the structure of thin film Nafion. Results show complex, layered interfacial structures, reduced water uptake, and reduced plasticity for films thinner than 60 nm due to confinement and substrate interactions [22,23]. However, there is limited experimental evidence directly connecting these structural properties and associated PEMFC CL fabrication/processing conditions to the thin film ionomer transport properties.

Detailed physical models provide insight into these losses and can help guide future PEMFC designs. Limiting processes in the CL occur over a wide range of length scales (ranging from nanometers to millimeter). Therefore, models typically implement effective cathode CL microstructures to capture relevant nm-scale phenomena within a continuum-scale model framework, using one of the three approaches: discrete-catalyst, thin-interface, and heterogeneous [14,24–26]. These types of models are mature in the PEMFC community, and a wide range of phenomena have been incorporated into them since the 1970s. The history of physically based PEMFC model development is well covered by a pair of reviews by Weber et al. [27,28]. Recent models have focused on phenomena such as the impact of microstructure [25,29,30], water management [30–33], Pt distribution [34,35], and species transport [32,36]. Broadly speaking, heterogeneous models offer the most detail and include effective microstructures such as the “core–shell” and “flooded agglomerate” models. As illustrated in Fig. 1, the core–shell microstructure takes as the computational domain a Pt-covered carbon core wrapped in an ionomer shell. The flooded agglomerate model considers a collection of Pt-covered carbon particles packed inside an ionomer shell that has been flooded with ionomer. Each approach has strengths and limitations. The core–shell model assumes that reactants and products only need to travel through a thin film of ionomer before reaching the gas phase. In reality, some of the Pt surface is“buried” within agglomerates, making reactant access more challenging. The flooded agglomerate model, on the other hand, underestimates the degree of Pt which is readily accessible to gas-phase reactants. The flooded agglomerate model requires artificially large agglomerates (hundreds of nanometers, encompassing tens of carbon support particles) in order to be compatible with continuum-scale finite differencing.

Actual microstructures, as observed in transmission electron microscope (TEM) images [37], lie somewhere between the two approaches. Agglomerates of roughly 3–5 carbon particles, flooded with ionomer, are typically observed. These agglomerates provide easy access to most of the Pt catalyst, but are typically highly non-spherical, with characteristic lengths ranging from roughly 50 to 300 nm. Such shapes are not readily represented in continuum-scale simulations, at present. Detailed computational fluid dynamics simulations can accurately represent high-fidelity microstructures, but scaling these up to the PEMFC stack level is a challenge. The effective geometries in Fig. 1, while not fully descriptive of the true CL microstructure, therefore provide reasonable computational performance at device-relevant scales.

In this work, we present pseudo-two-dimensional (2D) PEMFC simulations to explore limiting cathode CL phenomena at low Pt loading. Previously published structural data and ionic conductivity measurements are correlated with one another to develop structure–property relationships for the CL ionomer. These relationships capture transport losses related to variations in Pt loading, ionomer thickness, and CL microstructure, for both the core–shell and flooded agglomerate model frameworks. Results show that the core–shell geometry provides good agreement with experimental results for Pt loadings in the range of interest by the DOE. The flooded agglomerate models, on the other hand, require unrealistic ionomer transport properties and do not provide good agreement with data at low Pt loadings. Moreover, results demonstrate that microstructure models informed by TEM images and parameters based on physically meaningful structure–property relationships can capture complex transport and electrochemical kinetic phenomena in the PEMFC CL. With additional refinement, such modeling tools therefore provide a viable pathway to identify future PEMFC designs and operating conditions for low-cost, high-performance PEMFCs for EVs and other consumer applications.

## 2 Model Formulation

The model solves physically based conservation equations using a finite-volume approach to discretize the cathode GDL and CL. Within the GDL, one-dimensional porous flux expressions are used to compute evolution of the gas-phase species. For the CL, a pseudo-2D Newman-type model is used to track the PEMFC state variables as a function of the electrode depth (*y*) and the spherical solid particle radius (*r*). The mass density of oxygen absorbed in the Nafion ($\rho O2(Naf)$, $kgmNafion\u22123$) is tracked in the *y* and *r* directions. The surface site coverage of the Pt catalyst ($\theta k,surf$), when considered, is tracked in *y* and *r* (note that $\theta k,surf$ is not tracked in *r* for the core–shell model, since there is no radial distribution of Pt in the core–shell model; see below for details). The Nafion electric potential (Φ_{Naf}, V) and the gas-phase species mass densities ($\rho k,g$, $kgmgas\u22123$) are tracked only in the *y* dimension.

The model domain is illustrated in Fig. 1. Conservation equations are derived as time-based derivatives, solved numerically on a discretized series of finite volumes representing the model domain. The model equations capture the following physical processes: (i) gas-phase convective and diffusive transport, (ii) transport of ions and dissolved neutral species in the ionomer electrolyte, and (iii) reversible chemical and electrochemical reactions at ionomer-gas and ionomer-Pt interfaces. Similar model implementations are common in the literature [28], but what sets this model apart is the derivation and implementation of the ionomer structure–property relationships, as presented in Sec. 3.1. Evaluating the unsteady governing equations for these processes at a constant, user-defined current density and integrating over a long time span (100 s) until steady state provides the model output. Although the proton exchange membrane (PEM) is shown in Fig. 1, it is modeled here as a simple resistor with constant resistance. The anode is additionally not shown or explicitly modeled here, but its reversible electric potential is calculated and used as a reference so that the half-cell model can be directly compared to full-cell PEMFC experiments.

### 2.1 Model Assumptions.

The following assumptions are employed in the model:

Isothermal operation.

Water is in equilibrium between the Nafion and gas phases, and the gas remains at a fixed, steady relative humidity (RH).

Solid core–shell and/or flooded agglomerate structures are spherical in shape.

Uniform porosity, microstructure, and Pt distribution throughout the CL.

Electric potential is radially uniform within the spherical solid particles.

The carbon and Pt phases are isopotential through the CL depth.

Aside from neglecting fluctuations in water content via the second assumption, these assumptions are common throughout PEMFC modeling [32,38,39]. Although water transport and flooding are critical phenomena which limit cell performance with low Pt loading [31,40,41], they are neglected here so that the impact of structure–property relationships can be evaluated in the absence of confounding impacts from uncertain water dynamics. Incorporating more complex water dynamics and electrolyte phenomena will be the topic of future modeling work from our group.

### 2.2 Catalyst Layer Microstructure Modeling.

Two different structures are considered as representative domains within the CL: core–shell and flooded agglomerate. Each domain has been commonly used throughout PEMFC modeling literature [14,26]. Since the actual microstructure of a PEMFC CL is heterogeneous and non-uniform over very small length scales, these structures attempt to capture the average behavior within each discretized volume. The zoomed in section of the CL in Fig. 1 illustrates each domain. Note that the domains are not drawn to scale, but scale bars in the figure show the approximate relative sizes. A variety of values have been used with each of these geometries in previous literature, even beyond the scales provided here.

### 2.3 Conservation Equations

#### 2.3.1 Gas-Phase Mass and Species.

*k*is conserved by applying the expression

*A*

_{Naf}is the specific area of Nafion–gas interface per unit volume of CL. $\rho k,gas$,

*J*

_{k},

*W*

_{k}, and $s\u02d9k,Naf$ are the mass density ($kgkmgas\u22123$), mass flux ($kgkmCL\u22122s\u22121$), molecular weight ($kgmolk\u22121$), and net species molar production rate ($molkmNaf\u2212gas\u22122s\u22121$) due to reactions at the Nafion–gas phase interface for species

*k*, respectively. The mass flux in the gas phase captures transport via diffusion and convection as

*K*

_{g}is the gas-phase permeability (m

^{2}),

*μ*is the viscosity (Pa·s),

*P*is the pressure (Pa), $\tau fac,g$ is the tortuosity factor (estimated using a Bruggeman correlation [42] with an exponent of −0.5), and where

*D*

_{k}and

*Y*

_{k}are the diffusion coefficient (m

^{2}s

^{−1}) and mass fraction of species

*k*, respectively. Since the model is isothermal, the mass density of each species determines state of the gas phase, via the temperature, total phase density ($\rho gas=\u2211\rho k,gas$), and mass fractions ($Yk=\rho k,gas/\rho gas$).

#### 2.3.2 Electrolyte Mass and Species.

_{2}replacing the generalized species index

*k*. The new terms for $DO2(Naf)eff$, $s\u02d9O2,Pt$, and

*A*

_{Pt}represent the effective diffusion coefficient (m

^{2}s

^{−1}) of oxygen absorbed in Nafion, net molar production rate of oxygen ($molO2mint\u22122s\u22121$) from the reaction occurring at the Nafion–Pt interface, and specific area of Pt per unit CL volume, respectively. Since the gradients in oxygen occur radially due to the Pt reactions within the spherical domains, diffusion of absorbed oxygen in the

*y*direction is not considered. For any radial nodes in the electrolyte that do not contain Pt, $s\u02d9O2,Pt$ is zero. Likewise, the absorption reaction rate, $s\u02d9O2,Naf$, is nonzero only at the outermost radius where the Nafion–gas interface exists.

#### 2.3.3 Conservation of Charge.

*i*

_{ext}) equals the ionic current density into the membrane. Moreover, in each discretized CL volume, charge neutrality in the ionomer phase is preserved via the double layer current density

*i*

_{dl}($Amint\u22122$), which compensates for the non-zero sum of the divergence of the ionic current density plus the local Faradaic current density:

*i*

_{Far}($Amint\u22122$) is defined below in Sec. 2.4 and

*i*

_{io}is the ionic current density ($AmNafion\u22122$). The double layer current determines the double layer potential, ΔΦ

_{dl}, which is the difference between the electrolyte and electrode material electric potentials within the CL, ΔΦ

_{dl}= Φ

_{elyte,CL}− Φ

_{ed,CL}. This is modeled as capacitive in nature:

*Q*

_{dl}is the double layer charge ($Cmint\u22122$), which leads to the following for the rate of change of ΔΦ

_{dl}:

*C*

_{dl}and

*A*

_{dl}are the double layer capacitance ($Fmint\u22122$) and specific double layer area per unit volume of CL, respectively. The cathode potential is chosen as the reference potential (Φ

_{ca,CL}= 0 V), such that ΔΦ

_{dl}= Φ

_{elyte,CL}.

*i*

_{io}, is Ohmic in nature:

*σ*

_{io}is the ionic conductivity of the ionomer (S m

^{−1}) and $\epsilon Naf$ and $\tau fac,Naf$ are the Nafion volume fraction and tortuosity factor in the CL. In actual PEMFC CLs, some of the ionomer volume in $\epsilon Naf$ exists not as thin-film coatings or within flooded agglomerates, but rather as discrete, single-phase “globules” with higher ionic conductivity [43]. To compensate for this unmodeled heterogeneity, the tortuosity factor of the ionomer phase within the catalyst layer is set equal to unity.

### 2.4 Kinetics and Reactions.

^{−}(ca) is an electron in a cathode electronic conductor phase (Pt or C). Thermochemical reactions (e.g., reactions (8) and (10)) assume concentration-dependent reversible mass-action kinetics:

*ν*

_{k,j}, $\nu k,j\u2032$, and $\nu k,j\u2032\u2032$ are the net, forward, and reverse stoichiometric coefficients for species

*k*in reaction

*j*, and where $\nu k,j=\nu k,j\u2032\u2032\u2212\nu k,j\u2032$. [

*C*

_{ac,k}] is the “activity concentration” of species

*k*. This study assumes ideal thermodynamics (activity coefficients equal unity), such that the activity concentration equals the molar concentration, [

*C*

_{k}] (mol m

^{−3}for Nafion species, mol m

^{−2}for surface species), and

*k*

_{fwd,j}and

*k*

_{rev,j}are the forward and reverse reaction rate coefficients, whose ratio is dictated by thermodynamic reversibility. For the charge transfer reactions in either mechanism (e.g., reactions (9) and (11)), a concentration-dependent Butler–Volmer expression calculates

*i*

_{Far}:

*n*

_{elec}is the elementary electrical charge transferred by the charge transfer reaction,

*β*is the charge transfer symmetry factor,

*T*is the temperature (K), and $R\xaf$ is the universal gas constant (J mol

^{−1}K

^{−1}). The composition-independent portion of the exchange current density $i\u2218,j*$ for charge transfer reaction

*j*is calculated using Arrhenius parameters. The molar production rate for a species

*k*associated with a charge transfer reaction

*j*is

*j*reactions for a given interface.

### 2.5 Surface Coverages.

*k*, $\theta k$ evolve with time:

*k*. When using the one-step mechanism, the intermediate species that cover the Pt surface are not resolved, and Eq. (15) is neglected.

### 2.6 Membrane.

*R*

_{Naf}(Ω m

^{2}) can be used to describe the Ohmic losses across the membrane as

*R*

_{Naf}is on the order of what has been experimentally measured [44].

### 2.7 Anode.

*g*

_{rxn}is used to calculate the equilibrium potential

### 2.8 Cell Voltage.

_{ca}= 0, we get

### 2.9 Boundary Conditions.

*Y*

_{k}at the gas channel/GDL interface. Between the CL and membrane, two boundary conditions are applied. For all gas phase species

*k*, there is no flux in the

*y*direction into the membrane:

*r*

_{carbon}. For the flooded agglomerate model, the Pt–Nafion and carbon–Nafion interfaces exist at all radial locations except for the outermost shell. This moves the boundary condition of zero diffusive flux to

*r*= 0, which is treated as a symmetry plane.

## 3 Parameter Estimation

Based on TEM images and literature, a minimal number of parameters, all physically based, are used in the model. Briefly, the CL microstructure is defined by user inputs for the carbon particle radius, Pt loading, Pt radius, porosity, CL thickness, and Nafion shell thickness. In addition, the flooded agglomerate model requires additional inputs for the agglomerate radius and the volume fraction of carbon within the agglomerate. All other microstructure parameters are derived from these inputs. GDL and CL gas permeabilities *K*_{g} are taken from Ref. [29] and scaled linearly with the porosity. Input parameters that describe the core–shell and flooded agglomerate microstructures are summarized in Table 1. The parameters used in the model that are derived from this base set are included in Tables 2 and 3 for the core–shell and flooded agglomerate domains, respectively.

### 3.1 Structure–Property Relationships for Thin Film Catalyst Layer Ionomers.

Experiments show that transport properties for both protons and oxygen within the ionomer change as a function of the temperature, RH, thickness, and substrate interactions [21,45]. Specifically, variations in thickness and substrate type result in confinement effects and nano-structural rearrangements that impact water absorption and overall material structure [46,47]. Moreover, structural data from neutron reflectometry (NR) observe multi-layered Nafion structures at hydrophilic substrate interfaces [22,23,48]. These layers alternate between water-rich and water-poor compositions, with an oscillation that damps toward a uniform, thicker outer layer away from the substrate, beyond ca. 10 nm.

This model is unique in its use of structure–property relationships within the thin film CL ionomer, i.e., Nafion. These structure–property relationships are derived by combining three separate datasets:

*λ*:

*σ*

_{l}and

*f*

_{l}are the ionic conductivity and scaling factor for layer

*l*, respectively, and

*σ*

_{bulk}(

*λ*) is the reference ionic conductivity of a bulk Nafion membrane with the same

*λ*as the layer [49].

_{2}diffusion coefficients, experimentally measured $DO2$ values for thick films [50] are scaled by the local water content, relative to that in bulk membranes at similar RH, and are then multiplied by the transport scaling factor:

*λ*° is the water uptake of the reference membrane, and

*λ*

_{l}and

*f*

_{l}are the water uptake and scaling factor of layer

*l*, respectively.

The three structure–property methods used here make varying use of these scaling factors, as described as follows.

#### 3.1.1 Method 1: Layered.

The “layered” model hypothesizes that lamellae persist at every Nafion/solid interface. Effective diffusion coefficients for O_{2} absorbed in Nafion are calculated for each layer, as described in Eq. (24). Because the film is much thinner than the carbon particle radius, layer diffusion coefficients are summed using a linear series resistance to calculate an effective diffusion coefficient for the film. The ionic conductivities are taken directly from measurements on similar thickness films on silicon, as measured by McCreery et al. [21], interpolating when necessary.

#### 3.1.2 Method 2: Uniform.

Although layered interfacial structures have been observed on Si and Pt substrates using NR [22,23,51,52], similar definitive studies do not currently exist for carbon substrates. While it is possible that lamellae do not exist in PEMFC CLs, it is still known that confinement effects impact water absorption and transport through the Nafion electrolyte [22]. To account for this possible no-lamellae case, a “uniform film” method is implemented, which calculates uniform, isotropic transport properties for the thin Nafion films, using Eqs. (23) and (24), where *λ* is taken from NR measurements of similarly thick layers [22].

#### 3.1.3 Method 3: Mixed.

In all likelihood, if lamellae exist in the PEMFC CL, they are present to a limited extent. Layered structures have been observed at Nafion/Pt interfaces [51,52], consistent with the more hydrophilic Pt surface. The “mixed” model, then, proposes that lamellae are present at the Pt interface but that regions of exposed carbon are covered by “uniform” films. Because O_{2} only diffuses in regions near the Pt, $DO2,Naf,l$ for this model is calculated as in the “layered” approach. *σ*_{io}, meanwhile, is calculated as an average of the “uniform” and “layered” conductivities, with the two values weighted in parallel by the surface area fractions of carbon and Pt, respectively. The conductivity is therefore a function of the Pt loading. Higher loading correlates with higher conductivity, due to the higher water content of the interfacial layers.

### 3.2 Transport of Dissolved Nafion Species.

Figure 2 also demonstrates one additional degree-of-freedom for dissolved oxygen transport in the Nafion. For the core–shell model, a “capture angle” $\gamma $ is considered to represent lateral diffusion of the O_{2} within the Nafion film. Because the Pt does not cover the entire carbon surface, only a fraction of the Nafion coating will be active for O_{2} transport. If transport parallel to the film surface is negligible, the fraction of the Nafion area available for O_{2} transport at any given radius will be equal to the fraction of the carbon surface area covered by Pt. However, as “lateral” O_{2} transport becomes more prevalent, a greater fraction of the Nafion area becomes available. In reality, the transport field lines will be curved. As a first approximation, an effective capture angle $\gamma $ is used here to represent the phenomenon linearly, as shown in Fig. 2. Low $\gamma $ values represent the mostly radial diffusion, whereas larger $\gamma $ results in greater Nafion area available for O_{2} transport. Because O_{2} transport in the flooded agglomerate model is isotropic, with respect to the Pt surface orientation, this parameter is only implemented in the core–shell model.

## 4 Simulation Procedure

### 4.1 Numerical Integration.

Written as a system of ordinary differential equations in python,^{2} the model is solved numerically using the solve_ivp integrator available in the SciPy package [53]. A sufficiently large simulation time is used to ensure that steady-state conditions are reached. Running the model over a range of external currents between 0 and 2 A/cm^{2}, polarization curves are compared to experimental results with Pt loadings between 0.025 and 0.2 mg/cm^{2} from Ref. [54].

### 4.2 Cantera.

The open-source software library Cantera is used to manage the model’s thermodynamic, kinetic, and gas-phase transport parameter calculations [55]. Specifically, Cantera is used here to calculate mixture-averaged diffusion coefficients and viscosity for ideal gases, net production rates from the previously mentioned kinetic expressions, gas phase pressure and molecular weights. After setting the state of the Cantera phase objects using state variables, these properties can be calculated via simple and generalizable calls to Cantera library functions, simplifying the code structure significantly and increasing the repeatability of the model code. Cantera also makes it simple to change between the one-step and two-step charge transfer mechanisms, with no changes to the model code itself [56].

## 5 Accessing the Model

To advance open science practices related to accessibility, transparency, and repeatability, the model code is made freely available via GitHub [57]. We encourage any interested parties to download, use, adapt, or extend this model. Contributions to the model via GitHub are also welcome.

## 6 Results and Discussion

### 6.1 Pt-Loading Effects.

A goal of this model is to provide a physically meaningful explanation for fuel cell losses observed at low Pt loading. To accomplish this task, a set of five parameters including the kinetic rate coefficient for the surface adsorption step (Eq. (10)) *k*_{fwd,ads}, the concentration-independent portion of the exchange current density *i*_{o}*, membrane resistance *R*_{Naf}, Nafion O_{2} thermodynamics (characterized here by the free energy of the absorption reaction in Eq. (8), $\Delta gO2,abs$), and capture angle $\gamma $ are varied to find a set that, when held constant across all Pt loadings, provide a fit with relatively low error. The fitting routine utilizes the L-BFGS-B (i.e., the limited-memory Broyden–Fletcher–Goldfarb–Shanno algorithm with bounds) method of the ‘minimize’ package available within SciPy [53]. This algorithm minimizes the chi squared calculated between the model output and data, including variance, while allowing bounds for each parameter. The algorithm is in the family of quasi-Newton methods and uses an estimated hessian matrix to predict directions of improvement for each variable at each step. The optimized values from this fitting procedure are summarized in Tables 4 and 5. Figure 3 shows the best fits for both the core–shell and flooded agglomerate models. The core–shell model is shown using the “mixed” method for Nafion transport properties, as described in Sec. 3.1.3. Details for the properties used in the flooded agglomerate results are discussed below.

As observed in Fig. 3(a), the core–shell model provides a good fit across the range of Pt loadings, with slight discrepancies on the order of 10–20 mV for the lowest Pt loading. These discrepancies, in part, are likely due to the lack of modeling water transport. With the inclusion of liquid water, localized electrode flooding can limit transport at very low Pt loadings and high currents. Although the present study clearly identifies the importance of accurate structure–property relationships for the CL ionomer, incorporating water dynamics will likely improve the fit and change some of the fitted parameter values reported herein. This is the topic of forthcoming work from our group. Figure 3(b) demonstrates that the flooded agglomerate model is not capable of fitting the data. The fit captures trends for the highest two Pt loading values, but additional losses at low Pt are not captured. Moreover, obtaining the reported fit required abandoning the physically informed structure–property relationships developed herein. Using any of the three methods in Sec. 3.1 provides unrealistic limiting currents, related to oxygen transport limitations through the agglomerate. For the fit in Fig. 3, $DO2$ values of roughly 8 × 10^{−10} m^{2}/*s* were adopted from Ref. [25]. This value is two orders of magnitude larger than those generated from the structure–property relationship herein.

Comparing the results for the core–shell and flooded agglomerate models, it is apparent that physically realistic $DO2$ values are required to capture Pt loading effects. In the flooded agglomerate model, the high $DO2$ values compensate for the large agglomerate diameter. These large diameters, in turn, are required for valid application of continuum methods to the agglomerate interior. In other words, developing homogenized effective transport properties for strictly radial diffusion through a 150 nm agglomerate is not appropriate, when said agglomerate consists of only 2–5 carbon particles with diameter 50 nm. Such treatment requires larger agglomerates, which in turn must be compensated by artificially large $DO2$ values, for accurate limiting current predictions. However, these high $DO2$ values reduce the sensitivity of the model results to other important structure–property relationships, namely, the conductivity variations with Pt loading predicted by the “mixed” method. A higher affinity for liquid water near Pt surfaces (relative to carbon) correlates to higher ionic conductivities when Pt is more abundant. When $DO2$ is artificially high, the model is no longer sensitive to these variations in *σ*_{io}. Lastly, the core–shell geometry accommodates use of the effective capture angle, providing an additional relationship between the Pt loading and transport losses.

However, at the lowest Pt loading of 0.025 mg/cm^{2}, the core–shell results differ slightly from experimental data, demonstrating a general limitation of these effective microstructure approaches. In a real PEMFC CL, not all of the Pt is as easily accessible to oxygen as the core–shell geometry allows. At very low loading, Pt wedged between carbon cores likely has lower utilization, an effect that the core–shell model does not capture. A reasonable adaptation to hybridize the two microstructures presented here may provide a means to capture this effect. Use of a flooded agglomerate microstructure that neglects Nafion volume fraction and tortuosity effects and implements a transport capture angle for the outermost radial Pt-containing node could result in more realistic limiting currents, with the added benefit of capturing additional concentration losses at very low Pt loadings caused by limited access to oxygen. Alternately, higher fidelity simulations of a single, spatially resolved agglomerate can be used to derive effective transport parameters, which can then serve as inputs into reduced-order models with effective microstructure representations. The development of alternative modeling strategies is the subject of future work from our group.

Owing to its superior fit, the remainder of the results section will present results for only the core–shell model. This microstructure approach will provide a framework for analyzing other model phenomena, including structure–property approaches, Nafion thickness effects, and the impact of surface-mediated kinetic mechanisms.

### 6.2 Comparing Structure–Property Methods.

Of the three derived structure–property relationships, the “mixed” method provides the greatest sensitivity to Pt loading. Once operating conditions and shell thickness are chosen, the layered and uniform method transport properties are invariant with the Pt loading. The mixed method properties, on the other hand, are an average of those from the other two methods, weighted by the Pt loading. To examine the effects of these relationships, Fig. 4 shows predicted polarization curves for 0.05 mg/cm^{2} loading, using the layered and uniform structure–property methods. The effective ionic conductivities are 0.9 and 1.25 S/m for the uniform and layered methods, respectively, a difference that has little impact on the overall polarization behavior. Rather, the primary variable impacted by the structure–property method is the limiting current. The diffusion coefficient for oxygen in Nafion decreases by roughly 50% for the layered case, relative to the uniform model, from 6.3 × 10^{−12} to 3.2 × 10^{−12} m^{2}/s. This correlates with a limiting current of 2.8 A/cm^{2} in the layered model compared to 3.6 A/cm^{2} in the uniform model. The Fig. 4 inset demonstrates that this trend is consistent across Nafion shell thicknesses.

### 6.3 Shell Thickness Effects.

Figure 5 shows the polarization and power curves of the core–shell model at a Pt loading of 0.05 mg/cm^{2} with varying ionomer shell thicknesses. In this figure, as the Nafion shell thickness changes, the accompanying change in volume fraction is accommodated by changes in the gas-phase volume fraction. For the 25 nm case, the porosity is 0.1, but as the shell is reduced to 10 nm and 5 nm, the modeled CL porosity changes to approximately 0.5 and 0.6, respectively. This implementation ensures that the Pt surface area and carbon radius do not change so that any changes in polarization behavior are due purely to transport effects. It is clear from the figure that the 25 nm shell provides the highest losses, which are attributed to poor gas transport (due to low $\epsilon g$) and large O_{2} radial diffusive distances through the CL ionomer.

Another interesting result is that the 10 nm shell has a higher maximum power density than the 5 nm shell. This 13% increase in power density is despite lower porosity and the larger diffusive distances for the 10 nm case (which lead to a lower limiting current in the 10 nm case). The implementation of structure–property relationships in this model predicts higher ionic conductivity and oxygen diffusion coefficients for thicker shells due to greater water absorption. The model demonstrates that, within the range of realistic property co-variations, improved transport properties can outweigh correlated deleterious microstructural impacts. This result stands out from other models that assume transport properties which are invariant with Nafion thickness and Pt loading [25,26]. These other models generally use the shell thickness as a fitting parameter to adjust transport resistance. By making the shell thicker, they can capture additional resistance; however, in this model, thicker shells can result in less resistance due to improved transport properties.

To better understand predicted polarization trends as a function of CL Nafion thickness, Figs. 6 and 7 show profiles for the ionomer electric potential, Faradaic current distribution, and O_{2}(Naf) density as a function of CL depth. Results are presented for current density of 1.5 A/cm^{2} (i.e., where the performance of the 10 nm shell is noticeably better than that of the 5 nm case). Orientation along the *x*-axis is with the PEM interface at the left and the GDL interface at the right. To explain differences in ionic transport, Fig. 6 shows the electrolyte potential and Faradaic current distribution (fraction of total Faradaic current) throughout the CL depth. This plot demonstrates the impact of higher conductivities due to higher water absorption. The thicker 10 nm shell has higher ionic conductivity, such that utilization of protons is reasonably uniform throughout the CL depth, even near limiting currents. In the thinner 5 nm shell, we observe larger ionomer electric potential gradients due to lower *σ*_{io}, which in turn leads to a greater fraction of the Faradaic current near the electrolyte membrane. The larger Ohmic penalty is therefore associated with higher activation losses, in the 5 nm thick ionomer film case. The combination of losses provides an explanation as to why the thinner ionomer shell predicts lower performance, at this current density.

Looking at the oxygen depth profile in Fig. 7, we observe that the uneven Faradaic current of the 5 nm shell results in a similarly uneven distribution of O_{2}(Naf) concentration gradients throughout the CL depth. The oxygen density of the outermost shell is plotted for both film thicknesses to show that porosity effects do not significantly impact the gas-phase transport. Differences between the outermost and innermost shells represent the radial concentration gradients. For both the 5 and 10 nm film cases, we see that the Faradaic currents are transport-limited near the electrolyte membrane (O_{2}(Naf) densities are near zero at the Pt interface). However, even though the diffusion coefficient is higher for the 10 nm case, the O_{2}(Naf) gradients are more severe farther from the membrane, relative to the 5 nm case. This is due in part to the larger diffusive distance that oxygen travels to reach the reactive Pt sites in the 10 nm case but is also consistent with higher catalyst utilization farther from the membrane, in the 10 nm thickness case. Despite these larger radial diffusive gradients, which lead to a lower limiting current in the 10 nm case, this is more than compensated by the increase in ionic conductivity, leading to improved power density for a range of reasonable operating current densities.

### 6.4 Reaction Mechanism and Surface Chemistry.

All results to this point are based on the two-step reaction mechanism (Eqs. (10) and (11)). As shown in Fig. 3, this approach captures cell polarization phenomena in the low-current “activation region,” for both the core–shell and flooded agglomerate models over the range of Pt loadings. To explore the impact of detailed surface kinetics, a separate fit to the data was run with the one-step global mechanism (Eq. (9)). Fitting results show that Pt loading effects are not correctly predicted with this mechanism. Figure 8 shows a close-up of the kinetic region from polarization curves created with the core–shell geometry using the two reaction mechanisms. While the two-step mechanism, repeated in Fig. 8(a) for comparison’s sake, matches the data well within this region, Fig. 8(b) shows the one-step mechanism, at its best, results in a grouping of curves centered near the middle of the data. This grouping is caused by poorly predicted activation overpotentials resulting from simplified kinetics.

Equation (13) provides an explanation for this effect. In the Butler–Volmer formulation, the exchange current density is a function of the local activity concentrations [*C*_{ac,k}] (here, equal to the molar concentrations) of the products and reactants. In the one-step mechanism, the only such species whose activity concentration varies with current density is O_{2}(Naf). The proton and H_{2}O activity concentrations are constant, and so the exchange current density scales with [O_{2}(Naf)]^{0.5}. However, for the two-step mechanism, the intermediate species O(s) and Pt(s) participate directly in the reaction, with $\theta O(s)=1\u2212\theta Pt(s)$. Therefore, as the reactant O(s) is depleted at higher current densities, the product species Pt(s) concentration increases. The exchange current density for this mechanism scales with the ratio [O(s)]/[Pt(s)]. For fast, near-equilibrium kinetics, this ratio will be proportional to [*C*_{ac,k}], and the two mechanisms are roughly equivalent. However, for slow adsorption kinetics (Eq. (10)) associated with decreased Pt loading, the ratio [O(s)]/[Pt(s)] will drop faster than [O_{2}(Naf)]^{0.5} with increasing current. With decreasing Pt loading, the kinetic limitation becomes more severe, leading to the Pt loading sensitivity observed in Figs. 3 and 8(a), for the two-step mechanism. To state the obvious: while the global reaction is mathematically convenient and can accurately capture polarization in non-kinetic-limited performance, it is incumbent on the modeler to explore the implications of any kinetic mechanism assumptions. In a recent article, Dickinson and Hinds provide an insightful overview of various approaches to Butler–Volmer kinetics in PEMFCs and highlights the strengths and limitations of a number of common approaches [58].

## 7 Conclusions

To reduce the cost of PEMFCs while maintaining suitable performance, a quantitative and mechanistic understanding of the high observed overpotentials at low Pt loading is required. Two physically based models developed here identify the importance of both transport and kinetic phenomena in understanding and addressing these performance limitations. Physically informed structure–property relationships are hypothesized to describe the transport of ions and dissolved O_{2} through nano-thin Nafion sheets in the cathode CL. Functionally, these relationships only provide reasonable results for one of the two considered microstructures (the so-called “core–shell” model). With the use of (i) the core–shell model approach, (ii) transport properties based on experimentally observed Nafion structures which are sensitive to Pt loading, and (iii) a two-step, surface-mediated reaction mechanism, cell polarization behavior is matched over a range of low Pt loadings. This result provides a more detailed explanation for these losses compared to generalized resistances used in other models. Capturing Nafion transport properties based on the ionomer structure also provides insight into future design considerations. For example, a slight increase to the CL ionomer thickness predicts increased maximum power densities, due to improved ionic conductivity associated with greater water uptake in the CL Nafion. This work establishes that structure–property relationships based on ionomer thickness have non-negligible impact on the predicted performance and should be considered in future modeling efforts. Future development of new modeling frameworks for effective CL microstructures and incorporation of additional physical phenomena can lead to quantitatively predictive modeling tools to guide the design of advanced PEMFCs [34–36].

## Footnote

## Acknowledgment

This research was supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Materials Science and Engineering under Early Career Award #DE-SC0018019, Dr. Pappan Thiyagarajan, Program Manager.

## References

^{®}Membranes Measured by Neutron Reflectivity for PEM Fuel Cells