A numerical Assessment of Mitigation Strategies to Reduce Local Oxygen and Proton Transport Resistances in Polymer Electrolyte Fuel Cells

1. Introduction

Polymer electrolyte membrane fuel cells (PEMFCs) are one of the most promising alternatives to internal combustion engines due to their high efficiency, high power density, low noise and low pollutant emission (only water and heat are generated by electrochemical reactions) [1,2,3]. However, the high cost caused by the use of elevated Pt loadings at the cathode still remains as a critical issue for the widespread commercialization of PEMFCs either in portable, transport or stationary applications [4,5]. As a result, a large body of current research is devoted to mitigate potential losses in low Pt loading electrodes (see, e.g., [6,7,8,9] and references therein).

Operation at low Pt loading is limited by the sharp increase of the so-called local oxygen transport resistance around Pt sites, which originates a high concentration overpotential at the cathode [10,11]. The origin of the resistance is inherently caused by the reduction of active Pt surface area at low Pt loading, so that the local flux toward each catalyst particle is amplified to maintain a targeted current density (per unit of geometric or platform area). In other words, the lower the active surface area, the higher the oxygen flux that each catalyst particle must support [12]. Therefore, electrochemical activity and transport in the neighbourhood of every catalyst particle must be optimized to mitigate the limitation (i.e., high flux demand) caused by a reduction of Pt loading [7].

Previous work dealing with the design of CLs with improved performance at low Pt loading is reviewed below. In 2006-2012, Debe et al. [13,14] presented nanostructured thin film (NSTF) catalysts produced by 3M. Ultra-thin CLs (${\delta }_{cl}\simeq 0.3\mu \mathrm{m}$) were prepared with a Pt monolayer supported on crystalline organic whiskers without ionomer binder. Proton and electron conduction relayed on liquid water and Pt catalyst, respectively. Hih performance was reported a low Pt loading ($\simeq 0.1{\mathrm{mg}}_{\mathrm{Pt}}{\mathrm{cm}}^{-2}$) with superior stability in cyclic durability test compared to conventional Pt/C catalyst. The better performance at low Pt loading was ascribed to the absence of ionomer and the agglomerated Pt/whisker structure, which reduced Pt growth and coalescence due to Ostwald ripening. Yu et al. [15] (2015) examined the feasibility of reactive spray deposition technique for the preparation of high-performance CLs. They reported a better distribution of ionomer on active Pt sites and ionomer penetration into nanopores with a pore size of $1.7-10\mathrm{nm}$. As a result, the optimal ionomer-to-carbon weight ratio was decreased from around 0.65 for a conventional CL down to 0.3. The improved ionomer distribution and lower ionomer content led to higher performance at low Pt loading. Talukdar et al. [16] (2019) studied the effect of drying on the preparation of conventional ink-based CLs. They found that solvent removal via sublimation by means of freeze drying led to a 3.5-fold increase of the effective porosity, thereby enhancing performance at low Pt loading ($0.16{\mathrm{mg}}_{\mathrm{Pt}}{\mathrm{cm}}^{-2}$). Conde et al. [17] (2019) analyzed mass transport and water uptake properties of CLs prepared by electrospry. A significant decrease of the local oxygen transport resistance was reported compared to conventional CLs prepared by airbrushing ($0.8\mathrm{s}{\mathrm{cm}}^{-1}$ vs. $0.2\mathrm{s}{\mathrm{cm}}^{-1}$ at $0.025{\mathrm{mg}}_{\mathrm{Pt}}{\mathrm{cm}}^{-2}$, full humidification and $80{}^{\circ }\mathrm{C}$). Moreover, electrosprayed CLs showed a much higher water uptake than conventional CLs at micro-/meso-scale. Yoshino et al. [18] (2020) presented ionomer nanofiber scaffolding CL formed by the combination of electrospinning of Nafion nanofibers and electrospraying of a catalyst ink on nanofibers. The multiscale ionomer microstructure showed improved performance at low RH and reduced catalyst poisoning thanks to the decrease of the ionomer content in the catalyst ink. More recently, Cheng et al. [19] (2022) produced CLs with nanoporous ionomer using polyvinyl alcohol as a sacrificial pore-forming agent during CL preparation. Water-filled pores in ionomer allowed a dramatic reduction of the local oxygen transport resistance from $0.37\mathrm{s}{\mathrm{cm}}^{-1}$ down to $0.08\mathrm{s}{\mathrm{cm}}^{-1}$. Zhang et al. [20] (2022) developed porous ionomers by incorporating ionic covalent organic framework nanosheets into Nafion. Meso-porosity significantly enhanced oxygen permeation and mass activity, increasing performance at low Pt loading by a factor of 1.6 compared to conventional CLs.

The above literature review shows the key role that volume composition, morphology and ionomer porosity can play on CL performance. The aim of this work is to examine numerically the impact of CL microstructure on the local oxygen and ionic transport resistances around Pt sites. Predictions are compared against previous experimental observations to shed light on the various factors that affect performance at low Pt loading. The organization of this paper is as follows. In Section 2, the generation algorithm used for the virtual reconstruction of CL morphologies and the numerical model used for the calculation of local oxygen and ionic transport resistances are presented. The results are discussed in Section 3, which includes a calibration of the numerical model and a parametric analysis of volume composition, carbon/ionomer interaction, oxygen diffusivity in ionomer, and meso-porous ionomer. Finally, the conclusions are presented in Section 4.

2. Methods

2.1. Virtual reconstruction

CL microstructure was virtually reconstructed to differentiate among three phases, carbon, ionomer and free water, considering four ionomer morphologies: (i) heterogeneous, ($ii$) uniform, ($iii$) idealized vertically-aligned, and ($iv$) meso-porous ionomer distributions. The size of the representative domain was set to $300\times 300\times 300{\mathrm{nm}}^3$, significantly lower than average CL thickness, typically around ${\delta }_{cl}\sim 10\mu \mathrm{m}$. The voxel resolution was kept equal to $2\mathrm{nm}$ throughout the work. Figure 1 shows examples of the steps followed for the generation of heterogeneous CL microstructures, along with final modifications added to create meso-porous ionomer. Uniform and idealized morphologies were omitted in Figure 1 for brevity since they involve small changes of the algorithm used for the heterogeneous case. The main steps followed in the generation process were as follows:

• Carbon agglomeration. As shown in Figure 1(a), spherical carbon particles were randomly located in the domain with an uniform radius of $r_c=25\mathrm{nm}$, allowing overlapping between them. Particles were incorporated one by one into the domain until a prescribed carbon volume fraction, ${\epsilon }_c$, was reached. After each particle addition, only the largest connected component was maintained in the process, while isolated particles not connected to the main carbon structure were removed. At the end, the connectivity of the agglomerated carbon structure to the top and bottom surfaces of the domain was checked, and the generation process was repeated from the beginning if there was not a connected pathway across the domain (six-connected voxels criterion). Usually, no more than five iterations were needed to reach a connected structure at the lowest carbon volume fraction examined, ${\epsilon }_c=0.2$. The generation of the idealized carbon support composed of vertically-aligned cylinders was accomplished using a simplified algorithm. Carbon cylinders were placed with a uniform spacing in the material plane and the radius increased until reaching a prescribed carbon volume fraction.
• Ionomer addition. As shown in Figure 1(b), heterogeneous ionomer was created by randomly selecting points from the carbon agglomerate and introducing semi-spherical films around the structure. Ionomer films were incorporated one at a time by identifying the void voxels enclosed in an sphere centered at the selected carbon point with a prescribed ionomer radius, $r_i$. For each radius, ionomer films were sequentially added until no further variation of the ionomer volume fraction, ${\epsilon }_i$, was detected (below a established threshold). The whole process was completed when a prescribed porosity, $\epsilon$, was reached, gradually increasing $r_i$ by a factor of 1.2 from an exceedingly small value ($r_i=2\mathrm{nm}$). For uniform coating, the ionomer phase was simply identified using the Ecludiean distance transform, so that void voxels located at a distance below $r_i$ from the carbon pahse were identified as ionomer. As in the heterogeneous case, the ionomer radius was gradually increased by a factor of 1.2 from $r_i=2\mathrm{nm}$ until reaching the prescribed porosity, $\epsilon$. In all cases, connectivity was checked after every ionomer addition to remove isolated components not connected to the main carbon+ionomer structure.
• Free water addition. As shown in Figure 1(c), free water was added in a similar way to ionomer. However, random points were selected from either carbon, ionomer or water phases to identify void voxels to be converted into free water. The radius of water spheres, $r_w$, was sequentially increased by a factor of 1.2 from the last ionomer radius used in Step 2 until reaching a prescribed water saturation, s. In structures with uniform morphology, water was placed uniformly around ionomer by gradually increasing $r_w$ by a factor of 1.2. Isolated water blobs which were not connected to carbon or ionomer were removed.
• Additional features. Modifications of the morphology were incorporated to include specific features in the carbon, ionomer and water distributions. As shown in Figure 1(d), meso-porous ionomer was created by introducing water-filled spherical pores in the ionomer phase with a prescribed radius equal to the last ionomer radius used in Step 2. A total of 120 random points were selected from the ionomer phase.

2.2. Numerical model

Oxygen diffusion and ionic conduction were examined by direct numerical simulation of Laplace equation at pore scale (see, e.g., [21,22])

$$\nabla ·\left(\mathsf{\Gamma }\nabla \phi \right)=0$$

(1)

where $\mathsf{\Gamma }=D_{{\mathrm{O}}_2}^{\mathrm{eff}},{\sigma }_p^{\mathrm{eff}}$ is the property of interest, either mass diffusivity or ionic conductivity, and $\phi =C_{{\mathrm{O}}_2},{\varphi }_p$ is the variable of interest, either oxygen concentration or ionic potential.

The generated microstructures were imported into the CFD software ANSYS Fluent using a numerical mask to identify carbon, ionomer, free water and void regions in a hexahedral mesh created with the same voxel resolution of the image stacks. The various cell regions were then differentiated in the numerical mesh using the built-in capabilities of ANSYS fluent to separate cell regions and create the corresponding interfaces. Simulations were carried out on ionomer, water and void phases, considering a no-flux boundary condition at passive carbon interfaces were reaction does not take place. As shown in Figure 2, a Dirichlet boundary condition was set at the reactive carbon/ionomer interface, corresponding to limiting current density condition ($C_{{\mathrm{O}}_2}={\varphi }_p=0$). Additionally, the carbon/water interface was also assumed to be reactive ($C_{{\mathrm{O}}_2}={\varphi }_p=0$) in simulations with meso-porous ionomer. The system was excited by introducing a concentration or potential difference at the upper and lower surfaces ($C_{{\mathrm{O}}_2}={\varphi }_p=1$). Symmetry boundary conditions were set at sidewalls to mimic the representative location of the domain within larger carbon agglomerates.

The diffusivity and the ionic conductivity of ionomer and water were set constant. As for the void phase, diffusivity was corrected for Knudsen diffusion based on the average pore radius (determined by applying the Euclidean distance transform on image stacks). Ionic conductivity of the void phase was set to a negligible value due to the absence of conductive routes for protons. According to previous work, the values adopted in each region were as follows

$$D_{{\mathrm{O}}_2}^{\mathrm{eff}}=\left\{\begin{array}{cc}3\times {10}^{-11}{\mathrm{m}}^2{\mathrm{s}}^{-1}\hfill & \mathrm{ionomer}[23–26]\hfill \\ 4\times {10}^{-9}{\mathrm{m}}^2{\mathrm{s}}^{-1}\hfill & \mathrm{water}[10,27]\hfill \\ {\displaystyle \frac{D_{{\mathrm{O}}_2}^{\mathrm{bulk}}}{1+{\displaystyle \frac{D_{{\mathrm{O}}_2}^{\mathrm{bulk}}}{D_{{\mathrm{O}}_2}^{kn}}}}}\hfill & \mathrm{void}[10,12]\hfill \end{array}\right.$$

(2a)

$${\sigma }_p^{\mathrm{eff}}=\left\{\begin{array}{cc}1\mathrm{S}{\mathrm{m}}^{-1}\hfill & \mathrm{ionomer}[28–31]\hfill \\ 0.1\mathrm{S}{\mathrm{m}}^{-1}\hfill & \mathrm{water}[10,32,33]\hfill \\ {10}^{-3}\mathrm{S}{\mathrm{m}}^{-1}\hfill & \mathrm{void}\left(\mathrm{assumed}\right)\hfill \end{array}\right.$$

(2b)

Comparatively, the effective oxygen diffusivity in ionomer ($D_{{\mathrm{O}}_2,i}^{\mathrm{eff}}\sim {10}^{-7}{\mathrm{cm}}^2{\mathrm{s}}^{-1}$) is around five and two orders of magnitude lower than that in void space ($D_{{\mathrm{O}}_2,v}^{\mathrm{eff}}\sim {10}^{-2}{\mathrm{cm}}^2{\mathrm{s}}^{-1}$) and liquid water ($D_{{\mathrm{O}}_2,w}^{\mathrm{eff}}\sim {10}^{-5}{\mathrm{cm}}^2{\mathrm{s}}^{-1}$). So, bulk void space and liquid water lead to a negligible concentration drop compared to ionomer, which is where almost the entire oxygen transport resistance is originated. In contrast, the effective ionic conductivity of ionomer (${\sigma }_p^{\mathrm{eff}}\sim 1\mathrm{S}{\mathrm{m}}^{-1}$) is around one order of magnitude higher than that of liquid water (${\sigma }_p^{\mathrm{eff}}\sim 0.1\mathrm{S}{\mathrm{m}}^{-1}$). The increase of the effective ionic conductivity in ionomer (one order of magnitude) is not as pronounced as the decrease of the effective oxygen diffusivity (two orders of magnitude).

In the above expressions, $D_{{\mathrm{O}}_2}^{\mathrm{bulk}}$ and $D_{{\mathrm{O}}_2}^{kn}$ are the bulk and Knudsen diffusivities of oxygen, given by [10,12]

$$\begin{array}{cccc}& D_{{\mathrm{O}}_2}^{\mathrm{bulk}}\hfill & & =2.65\times {10}^{-5}{\left({\displaystyle \frac{T}{333.15}}\right)}^{1.5}\hfill \end{array}$$

(3a)

$$\begin{array}{cccc}& D_{{\mathrm{O}}_2}^{kn}\hfill & & ={\displaystyle \frac{r_v}{3}}\sqrt{{\displaystyle \frac{8RT}{\pi M_{{\mathrm{O}}_2}}}}\hfill \end{array}$$

(3b)

where $r_v$ is the average pore radius of void space, T is the temperature (fixed to $80{}^{\circ }\mathrm{C}$), $M_{{\mathrm{O}}_2}$ is the molecular mass of oxygen, and R is the universal gas constant.

The local mass and ionic transport resistances per unit of geometric (or platform area of the CL), $R_{{\mathrm{O}}_2}^{\mathrm{local}}$ and $R_p^{\mathrm{local}}$, are obtained from a flux balance [12]

$$\begin{array}{cccc}& N_{{\mathrm{O}}_2,\mathrm{Pt}}^{\mathrm{local}}{r}_f={\displaystyle \frac{\Delta C_{{\mathrm{O}}_2}}{R_{{\mathrm{O}}_2}^{\mathrm{local}}}}\hfill & & \Rightarrow R_{{\mathrm{O}}_2}^{\mathrm{local}}={\displaystyle \frac{\Delta C_{{\mathrm{O}}_2}}{N_{{\mathrm{O}}_2,\mathrm{Pt}}^{\mathrm{local}}}}{\displaystyle \frac{1}{r_f}}\hfill \end{array}$$

(4a)

$$\begin{array}{cccc}& N_{p,\mathrm{Pt}}^{\mathrm{local}}{r}_f={\displaystyle \frac{\Delta {\varphi }_p}{R_p^{\mathrm{local}}}}\hfill & & \Rightarrow R_p^{\mathrm{local}}={\displaystyle \frac{\Delta {\varphi }_p}{N_{p,\mathrm{Pt}}^{\mathrm{local}}}}{\displaystyle \frac{1}{r_f}}\hfill \end{array}$$

(4b)

where $r_f=A_{\mathrm{Pt}}/A_{\mathrm{geo}}$ is the roughness factor, defined as the ratio between the active Pt surface area (ECSA) and the geometric area of the CL, $N_{{\mathrm{O}}_2,\mathrm{Pt}}^{\mathrm{local}}$ and $N_{p,\mathrm{Pt}}^{\mathrm{local}}$ are the oxygen and proton flux at Pt surface, and $\Delta C_{{\mathrm{O}}_2}$ and $\Delta {\varphi }_p$ are the concentration and ionic potential drop from the bulk space toward Pt surface, respectively. $\Delta C_{{\mathrm{O}}_2}$ was determined as the average concentration in the void space minus the oxygen concentration at reaction sites ($\Delta C_{{\mathrm{O}}_2}=C_{{\mathrm{O}}_2}^{\mathrm{bulk},\mathrm{avg}}$), while $\Delta {\varphi }_p$ was directly taken as the prescribed potential drop ($\Delta {\varphi }_p=1$). According to Eqs. (4a)–(4b), the local resistance per unit of geometric area is inversely proportional to the roughness factor ($R^{\mathrm{local}}\propto r_f^{-1}$), dramatically increasing at low Pt loading when $A_{\mathrm{Pt}}\to 0$ [11,12,34].

Pt was assumed to be well dispersed over Pt surface, so not explicit description of Pt nano-particle shape was made ($A_c\simeq A_{\mathrm{Pt}}$). Consequently, the roughness factor is meaningless, and the variable of interest calculated in the simulations was the local resistance per unit of Pt (or carbon) surface area $R_{\mathrm{Pt}}^{\mathrm{local}}$

$$\begin{array}{cccc}& R_{{\mathrm{O}}_2,\mathrm{Pt}}^{\mathrm{local}}\hfill & & =R_{{\mathrm{O}}_2}^{\mathrm{local}}{r}_f={\displaystyle \frac{\Delta C_{{\mathrm{O}}_2}}{N_{{\mathrm{O}}_2,\mathrm{Pt}}^{\mathrm{local}}}}\hfill \end{array}$$

(5a)

$$\begin{array}{cccc}& R_{p,\mathrm{Pt}}^{\mathrm{local}}\hfill & & =R_p^{\mathrm{local}}{r}_f={\displaystyle \frac{\Delta {\varphi }_p}{N_{p,Pt}^{\mathrm{local}}}}\hfill \end{array}$$

(5b)

$R_{\mathrm{Pt}}^{\mathrm{local}}$ controls the incremental slope of $R^{\mathrm{local}}$ when the Pt loading (and $A_{\mathrm{Pt}}$) are decreased, removing the expected effect of the increase of $N_{\mathrm{Pt}}^{\mathrm{local}}$ with respect to $N^{\mathrm{local}}$ caused by area (i.e., $N_{\mathrm{Pt}}^{\mathrm{local}}=r_f{N}^{\mathrm{local}}$). The only concern is that the enlargement of transport pathways created by the discrete nature of Pt nano-particles at a given Pt loading is ignored [35]. Nevertheless, the aim of this work is to examine the overall effect of ionomer morphology, so the contribution of the discrete Pt shape prevails as a secondary factor.

3. Discussion of results

The discussion of results is divided into five sections. In Section 3.1, the model is calibrated in terms of local oxygen transport resistance and effective diffusivity of oxygen in ionomer. The results of a parametric analysis are then presented, which includes a study of volume composition (Section 3.2), carbon/ionomer interaction (Section 3.3), ionomer diffusivity (Section 3.4), and meso-porous ionomer (Section 3.5). The local oxygen and ionic transport resistances computed in the parametric analysis were averaged among five sample realizations for every condition examined. Water saturation was fixed to $s=0.3$, while porosity and carbon volume fraction were varied depending on the case under study.

3.1. Calibration

Before proceeding further, the predictive capabilities of the model were compared against data reported in the literature. The variables considered for the model calibration were the local oxygen transport resistance, $R_{{\mathrm{O}}_2}{r}_f$, and the effective diffusivity of oxygen in ionomer, $D_{{\mathrm{O}}_2,i}^{\mathrm{eff}}$, due to the uncertainty in mass transport results and the availability of previous experimental data. Figure 3(a) shows the experimental data of $R_{{\mathrm{O}}_2}{r}_f$ reported by different literature sources, along with the numerical results computed with the baseline $D_{{\mathrm{O}}_2,i}^{\mathrm{eff}}$ used in the simulation campaign ($D_{{\mathrm{O}}_2,i}^{\mathrm{eff}}=3\times {10}^{-7}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1}$). A partially saturated CL ($s=0.3$) with a conventional composition (${\epsilon }_c=0.3$, $\epsilon =0.4$) was used in the simulations. A large dispersion is found among previous data, with $R_{{\mathrm{O}}_2}{r}_f$ ranging between $5\mathrm{s}{\mathrm{cm}}_{\mathrm{Pt}}^{-1}$ and $30\mathrm{s}{\mathrm{cm}}_{\mathrm{Pt}}^{-1}$ (a factor of six). These variations can be mainly ascribed to microstructural differences that arise from volume compositions and manufacturing techniques. The numerical results computed in sixteen sample realizations vary stochastically around $12\mathrm{s}{\mathrm{cm}}_{\mathrm{Pt}}^{-1}$, which is close to the mean experimental value ($\approx 15\mathrm{s}{\mathrm{cm}}_{\mathrm{Pt}}^{-1}$). Hence, the numerical results can be considered representative of the behavior of a conventional CL when the oxygen diffusivity is fixed to $D_{{\mathrm{O}}_2,i}^{\mathrm{eff}}\approx 3\times {10}^{-7}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1}$. As shown in Figure 3(b), the assumed diffusivity is in line with previous data adopted in mesoscopic simulations by Mu et al. [23,24], who varied the oxygen diffusivity between $2.5\times {10}^{-7}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1}$ and $1\times {10}^{-6}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1}$. This range agrees with previous experimental data reported for the oxygen diffusivity in bulk Nafion membranes, even though there is a large variability among authors (three orders of magnitude) [25]. The large fluctuation of the diffusivity in bulk Nafion can be ascribed to conditioning, substrate interaction and confinement of membranes [36,37,38]. The oxygen diffusivity considered here is a good approximation, which leads to realistic values of the local oxygen transport resistance.

3.2. Volume composition

Figure 4 shows the results of the analysis of the volume composition for conventional CLs with a heterogeneous ionomer distribution. The variation of the local oxygen and ionic transport resistances is examined as a function of porosity for three carbon volume fractions. As shown in Figure 4(a), there is a strong dependency of $R_{{\mathrm{O}}_2}{r}_f$ with porosity, which decreases by a factor of two when the porosity is increased from $\epsilon \simeq 0.1$ to $\epsilon \simeq 0.7$. That is, the ionomer volume fraction is reduced from ${\epsilon }_i\simeq 0.7$ to ${\epsilon }_i\simeq 0.1$ since $\epsilon +{\epsilon }_i+{\epsilon }_c=1$. For a given porosity, $R_{{\mathrm{O}}_2}{r}_f$ is also reduced with carbon volume fraction due to the decrease of ionomer volume fraction. These results agree with previous works, which showed that the use of moderate ionomer-to-carbon weight ratios ($\mathrm{I}/\mathrm{C}\lesssim 0.65$) is beneficial to enhance oxygen transport at low Pt loading due to a decrease of average ionomer thickness [12,15,34,43,46]. Notice also the non-linear dependency of the local oxygen transport resistance with porosity. For high ionomer volume fractions (${\epsilon }_i\simeq 0.5-0.7$), $R_{{\mathrm{O}}_2}{r}_f$ remains large and rather constant around $R_{{\mathrm{O}}_2}{r}_f\gtrsim 15\mathrm{s}{\mathrm{cm}}_{\mathrm{Pt}}^{-1}$. In contrast, at lower ionomer volume fractions (${\epsilon }_i\simeq 0.1-0.5$), $R_{{\mathrm{O}}_2}{r}_f$ gradually decreases until it settles down around $R_{{\mathrm{O}}_2}{r}_f\simeq 10\mathrm{s}{\mathrm{cm}}_{\mathrm{Pt}}^{-1}$ for CLs with low ionomer volume fractions and high porosities (${\epsilon }_i\simeq 0.1-0.3$, $\epsilon \simeq 0.5-0.7$). Such non-linear behavior highlights the need to design high porosity CLs to mitigate the detrimental effect of local oxygen transport resistance at low Pt loading [47]. Experimentally, a strong reduction of $R_{{\mathrm{O}}_2}{r}_f$ has been previously reported for high porosity electrosprayed CLs [17,48]. A low ionomer content minimizes the agglomeration of ionomer films, avoiding the formation of locally dense ionomer regions that prevent a proper distribution of oxygen throughout the ionomer surface. CL design with low ionomer fraction and high porosity must be accompanied of a moderate increase of the local ionic transport resistance. As shown in Figure 4(b), $R_p^{\mathrm{local}}{r}_f$ remains rather constant around $1.5-2\mathrm{m}\mathsf{\Omega }{\mathrm{cm}}_{\mathrm{Pt}}^{-2}$ for $\epsilon \lesssim 0.5$ (${\epsilon }_i\gtrsim 0.2$) but significantly increases beyond $3\mathrm{m}\mathsf{\Omega }{\mathrm{cm}}_{\mathrm{Pt}}^{-2}$ for $\epsilon \simeq 0.7$ (${\epsilon }_i\simeq 0.1$) when ionomer interconnection approaches the percolation threshold [49]. Operation close to the ionic percolation threshold of ionomer must be ensured by the incorporation of supporting routes for proton transport, created, for example, by an enhanced water uptake at micro- and meso-scale, as it is the case of optimized electrosprayed CLs [17,48]. Alternatively, a delicate arrangement of ionomer might be necessary, as discussed in the next section.

3.3. Carbon/ionomer interaction

Figure 5 shows the results of the carbon/ionomer interaction for three ionomer morphologies, heterogeneous, uniform and idealized vertically-aligned distributions, as a function of porosity. A representative carbon volume fraction is considered, ${\epsilon }_c=0.3$. As shown in Figure 5(a), the local oxygen transport resistance of the heterogeneous and uniform ionomer distributions are similar at low-to-middle porosities ($\epsilon \simeq 0.1-0.4$). The low sensitivity of $R_{{\mathrm{O}}_2}{r}_f$ to ionomer heterogeneity is explained by the presence of both thick and thin ionomer regions in the heterogeneous samples, which lead on average to a similar oxygen resistance to that of uniform samples [50,51,52]. This situation is further illustrated in Figure 6(a), which shows the histograms of the local oxygen flux at active sites of the examined samples at $\epsilon =0.4$. For heterogeneous ionomer, the local oxygen flux distribution varies roughly linear due to the presence of a continuous distribution of regions with high and low oxygen resistances. In contrast, for uniform ionomer, the local oxygen flux distribution is more concentrated toward intermediate values that are characteristic of the mean ionomer thickness (lower dispersion). The difference between heterogeneous and uniform ionomer distributions is enlarged at high porosity ($\epsilon \simeq 0.5-0.6$), where the uniform samples reach a lower oxygen resistance than the heterogeneous samples (a factor of two, $R_{{\mathrm{O}}_2}{r}_f\simeq 5\mathrm{s}{\mathrm{cm}}_{\mathrm{Pt}}^{-1}$ vs. $R_{{\mathrm{O}}_2}{r}_f\simeq 10\mathrm{s}{\mathrm{cm}}_{\mathrm{Pt}}^{-1}$). The reduction of the oxygen resistance at high porosity in the uniform samples is caused by a strong suppression of ionomer agglomeration when the ionomer volume fraction is low (${\epsilon }_i\simeq 0.1-0.2$) [43]. The optimal arrangement of uniform ionomer virtually removes coalescence of ionomer films, facilitating oxygen distribution between ionomer films. Similarly, the local ionic transport resistance of the uniform samples is especially reduced compared to that of the heterogeneous samples at high porosity (a factor of two, $R_p{r}_f\simeq 1.5\mathrm{m}\mathsf{\Omega }{\mathrm{cm}}_{\mathrm{Pt}}^{-2}$ vs. $R_p{r}_f\simeq 3\mathrm{m}\mathsf{\Omega }{\mathrm{cm}}_{\mathrm{Pt}}^{-2}$). The superior ionic conductivity of uniform samples arises from the better arrangement of ionomer and the suppression of ineffective agglomerated ionomer regions, i.e., a tortuosity decrease [53]. Overall, the use of high porosity CLs with a uniform as possible ionomer distribution can be an effective way to improve performance at low Pt loading due to facilitated oxygen transport and good proton conduction [46,54,55].

The importance of ionomer morphology is exacerbated in the samples with an idealized vertically-aligned distribution. As shown in Figure 5(a), the local oxygen transport resistance of the idealized samples is low in the full porosity range, slightly decreasing from $5\mathrm{s}{\mathrm{cm}}_{\mathrm{Pt}}^{-1}$ to $2\mathrm{s}{\mathrm{cm}}_{\mathrm{Pt}}^{-1}$ when the porosity is increased between $0.1$ and $0.6$. Furthermore, as shown in Figure 5(b), the local ionic transport resistance is significantly reduced compared to the heterogeneous and uniform distributions, remaining between $0.6\mathrm{m}\mathsf{\Omega }{\mathrm{cm}}_{\mathrm{Pt}}^{-2}-1.1\mathrm{m}\mathsf{\Omega }{\mathrm{cm}}_{\mathrm{Pt}}^{-2}$ for $\epsilon =0.1-0.6$. The extremely high performance of the idealized ionomer distribution arises from a minimization of ionomer agglomeration and an optimal arrangement of proton transport pathways (i.e., tortuosity minimization) [10]. Notice that the oxygen resistance of the idealized and the uniform distributions are similar at high porosity, reflecting the superior oxygen distribution in both cases. High performance at low Pt loading with vertically-aligned microstructures have been previously demonstrated in several experimental works [56,57,58], where current densities as high as $2.6\mathrm{A}{\mathrm{cm}}^{-2}$ were reached at $0.6\mathrm{V}$ with $0.1{\mathrm{mg}}_{\mathrm{Pt}}{\mathrm{cm}}^{-2}$ and air supply. As a final remark, it is worth noting that the results of the idealized morphology do not vary among sample realizations, unlike the results calculated with random carbon particles. This fact highlights the importance of using reproducible manufacturing techniques and homogenized catalyst supports to reduce maldistribution of transport resistances, which can affect local performance and degradation rates [59,60].

3.4. Ionomer diffusivity

The results of the analysis of the effective oxygen diffusivity in ionomer are shown in Figure 7. The local oxygen and ionic transport resistances are plotted as a function of porosity for the heterogeneous ionomer distribution (${\epsilon }_c=0.3$). As expected, the local oxygen transport resistance decreases almost linearly with an increase of the oxygen diffusivity (see Figure 7(a)). Indeed, small deviations from linearity are caused by stochastic fluctuations that arise from random microstructures. As shown in Figure 7(b), the local ionic transport resistance remains invariant with respect to oxygen diffusivity, reaching similar values to those presented in Figure 4(b) (except for stochastic variations among samples). The use of modified ionomers with enhanced permeability and good conductivity offers a viable route to improve performance at low Pt loading. Among available options, blended ionomers can be a satisfactory solution to produce tailored designs with balanced mass and ionic transport properties [61,62]. Nevertheless, it is worth noting that achieving a remarkable high performance may only be possible if enhanced diffusion is accompanied by a proper CL design in terms of high porosity and ionomer morphology. The increase of oxygen diffusivity shown in previous work is typically of order unity, so a disruptive increase of oxygen diffusivity may not be possible by solely altering the chemical ionomer structure [44,61,62,63,64]. Local oxygen transport can be more effectively increased by incorporating meso-porosity into ionomer, an aspect that is discussed in the next section.

3.5. Porous ionomer

The results of the uniform meso-porous ionomer are compared to those of the uniform ionomer as a function of porosity in Figure 8 (${\epsilon }_c=0.3$). As shown in Figure 8(a), the local oxygen transport resistance of the meso-porous ionomer is sharply decreased, reaching values even one order of magnitude lower than the uniform ionomer ($R_{{\mathrm{O}}_2}{r}_f\simeq 0.1-1\mathrm{s}{\mathrm{cm}}_{\mathrm{Pt}}^{-1}$ vs. $R_{{\mathrm{O}}_2}{r}_f\simeq 5-10\mathrm{s}{\mathrm{cm}}_{\mathrm{Pt}}^{-1}$). The strong enhancement of oxygen transport is caused by the much higher diffusivity of oxygen in liquid water ($D_{{\mathrm{O}}_2,w}^{\mathrm{eff}}\sim {10}^{-5}{\mathrm{cm}}^2{\mathrm{s}}^{-1}$) compared to ionomer ($D_{{\mathrm{O}}_2,i}^{\mathrm{eff}}\sim {10}^{-7}{\mathrm{cm}}^2{\mathrm{s}}^{-1}$), two orders of magnitude higher –see histograms in Figure 6(b). This result clearly evidences that oxygen resistance can be greatly suppressed by enhancing water uptake within ionomer and allowing oxygen transport directly through liquid water (provided that water regions are electrochemically activated). For instance, the increase of the oxygen resistance of the meso-porous ionomer with CL porosity is explained by a reduction of the carbon/water interface since ionomer pore radius was varied according to the film thickness (a decreasing function of CL porosity). As shown in Figure 8(b), the local ionic transport resistance shows an opposite behaviour to that of the oxygen resistance, being higher in the meso-porous ionomer due to the worse transport of protons in liquid water (a factor around two to three, $R_p^{\mathrm{local}}{r}_f\approx 1.5-4.5\mathrm{m}\mathsf{\Omega }{\mathrm{cm}}_{\mathrm{Pt}}^{-2}$ vs. $R_p^{\mathrm{local}}{r}_f\approx 1-2\mathrm{m}\mathsf{\Omega }{\mathrm{cm}}_{\mathrm{Pt}}^{-2}$). The ionic resistance of the meso-porous ionomer decreases with porosity due to the presence of smaller water regions at high porosity (i.e., narrower ionomer film thickness). Comparatively, the increase of the ionic transport resistance of the meso-porous ionomer is significantly lower than the reduction of the oxygen transport resistance, so reaching high performance with meso-porous ionomer is certainly possible. Previous experimental works that have shown improved performance at low Pt loading by exploiting the use of liquid water include: (i) nano-structured thin film electrodes [13,14], and ($ii$) engineered porous ionomers [19,20,65,66]. In the latter group, it may also be included some blended ionomers with a significant bulky structure [61]. In agreement with the present results, a reduction of the local oxygen transport resistance down to one order of magnitude ($R_{{\mathrm{O}}_2}=0.08\mathrm{s}{\mathrm{cm}}^{-1}$ vs. $R_{{\mathrm{O}}_2}=0.37\mathrm{s}{\mathrm{cm}}^{-1}$) with a moderate ionic resistance in the polarization curve was recently reported by Cheng et al. [19] using a masked nanoporous ionomer. In addition, the low oxygen transport resistance of electrosprayed CLs at low Pt loading ($R_{{\mathrm{O}}_2}^{\mathrm{local}}\simeq 0.17\mathrm{s}{\mathrm{cm}}^{-1}$) can be explained by the combination of high porosity, more uniform ionomer distribution and enhanced water uptake at micro-/meso-scale [17,48].

4. Conclusions

The local and ionic transport resistances of catalyst layers (CLs) used in polymer electrolyte membrane fuel cells (PEMFCs) have been analyzed by direct numerical simulation on virtually reconstructed microstructures. The generated samples accounted for carbon, ionomer and liquid water, neglecting the discrete geometry of Pt nano-particles, to study the effect of volume composition, carbon/ionomer interaction, ionomer diffusivity and ionomer meso-porosity. To this end, four ionomer morphologies were considered: (i) heterogeneous, ($ii$) uniform, ($iii$) idealized vertically-aligned, and ($iv$) meso-porous distributions. The numerical results have shown that the local oxygen transport resistance can be reduced through the design of high porosity CLs with low ionomer agglomeration. Increasing the uniformity and reducing the tortuosity of ionomer distribution avoids excessive coalescence of ionomer films, while enhancing proton transport. Oxygen transport can be further improved through the use of blended ionomers with high oxygen permeability and good ionic conductivity, as well as meso-porous ionomers. Allowing oxygen diffusion through liquid water near active Pt sites is crucial to strongly enhance oxygen transport toward Pt surface. Indeed, the high diffusivity of oxygen in liquid water compared to ionomer can decrease the oxygen resistance by one order of magnitude. The above guidelines, high porosity, uniform and low tortuosity ionomer and oxygen transport through activated liquid water, are key aspects to mitigate losses at low Pt loading in next-generation PEMFCs.

Several aspects warrant future work. The effect of the discrete geometry of Pt nano-particles and porous carbon support must be examined. In addition, a numerical analysis accounting for electrochemical reactions must be performed, with special emphasis on reaction kinetics in liquid water.