Numerical Study of the Impact of Two-Phase Flow in the Anode Catalyst Layer on the Performance of Proton Exchange Membrane Water Electrolysers

The governing equations are:

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot \left({\sigma }_{m}^{{eff}}{\rm{\nabla }}{\phi }_{m}\right)={S}_{{H}^{+}}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot \left({\sigma }_{s}^{{eff}}{\rm{\nabla }}{\phi }_{s}\right)={S}_{{e}^{-}}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot \left(\displaystyle \frac{{\rho }_{g}{k}_{g}}{{\mu }_{g}}{\rm{\nabla }}{p}_{g}\right)={S}_{{gas}}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot \left(\displaystyle \frac{{\rho }_{l}{k}_{l}}{{\mu }_{l}}{\rm{\nabla }}{p}_{l}\right)={S}_{{liquid}}\end{eqnarray} \tag{ 4 }$$

where the source terms S_i for each equation are based on the current produced in each catalyst layer, as shown in Table I and is discussed below. The effective layer electronic conductivity is denoted ${\sigma }_{s}^{{eff}}$, while the protonic is ${\sigma }_{m}^{{eff}}$. k_i is the permeability, μ_i is viscosity and ρ_i is the density of phase i, which is either water or oxygen. Note that Zhou et al. ⁵⁵ considered the evaporation and condensation of water as additional source terms for the liquid and gas pressure equations, however it is assumed in this work that the oxygen bubbles formed in the anode are fully humidified. Many of the assumptions made in reference ⁵⁵ also apply in this work, namely: liquid transport is dominated by surface tension forces, there exists a representative elementary volume in each layer that contains homogeneous hydrophilic and hydrophobic pore networks, gas and liquid water are incompressible and the membrane is impermeable to water and gas. Note that because of the membranes impermeability, liquid water cannot exist on the cathode side. It is assumed that the cathode electrode does not contribute significantly to the cell performance and so two-phase flow is not studied there. The permeability of the liquid, k_l , and gas phase, k_g , is computed using the PSD model as is discussed below. The viscosity, μ_g , and density, ρ_g , of the gas phase are computed assuming a mixture of oxygen that is fully saturated with water vapor.

Table I. Source terms for the governing equations.

Source term	Anode	Cathode
${S}_{{H}^{+}}$	jOER	−jHER
${S}_{{e}^{-}}$	−jOER	jHER
Sgas	${j}_{{OER}}{M}_{{O}_{2}}/4F$	$-{j}_{{HER}}{M}_{{H}_{2}}/2F$
Sliquid	$-{j}_{{OER}}{M}_{{H}_{2}O}/2F$	—

The kinetics of the OER and hydrogen evolution reaction (HER) are typically described by the Butler-Volmer or Tafel equation in the modelling literature, ⁶¹ with the exception of Oliveira et al. ^62,63 who developed a multi-scale cell model with multi-step kinetics for the OER. The nature of the OER is not well understood, with numerous mechanisms proposed that have not been experimentally validated. ²¹ The primary concern for this work is a model that can capture any change in the Tafel slope for the OER, which is occasionally seen experimentally. ^23,24 Ma et al. ²⁴ proposed a simplified reaction mechanism using the same methodology proposed by Wang et al. ⁶⁴ for the ORR. The model assumes a single intermediate species with two intermediate reactions, and was used by the authors to capture the kinetics of the OER on Ir/Ru oxide core-shell catalysts. This model is used in this work to capture the reaction kinetics as a Tafel model was not able to provide an accurate fit to the experimentally obtained polarization curves as is discussed below. The model assumes the following intermediate steps:

$$\begin{eqnarray}&&1/2\ {H}_{2}O\rightleftarrows 1/2\ {O}_{{ad}}+{H}^{+}+{e}^{-}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&1/2\ {H}_{2}O+1/2\ {O}_{{ad}}\rightleftarrows 1/2\ {O}_{2}+{H}^{+}+{e}^{-}\end{eqnarray} \tag{ 6 }$$

from which the current density can be computed using:

$$\begin{eqnarray}&&{j}_{{OER}}={s}_{L}{A}_{v}{j}^{* }\left({g}_{+1}(1-{\theta }_{O})-{g}_{-1}{\theta }_{O}\right)\end{eqnarray} \tag{ 7 }$$

where the active area is given by A_v and the reference prefactor, j^*, is used to scale the reaction. In practice, it was found that this model was also not able to capture the kinetics of the OER for the TKK catalyst, as very low Tafel slopes were required, necessitating the use of transfer coefficients greater than one and therefore forward and backward transfer coefficients that are decoupled (i.e. the backward is not one minus the forward). Then, the coverage of the assumed intermediate species, θ_O , is given as:

$$\begin{eqnarray}&&\theta =\displaystyle \frac{{g}_{+1}+{g}_{-2}}{{g}_{+1}+{g}_{-1}+{g}_{+2}+{g}_{-2}}\end{eqnarray} \tag{ 8 }$$

and the g_i terms are given as:

$$\begin{eqnarray}&&{g}_{+1}=\exp \left(-\displaystyle \frac{{\rm{\Delta }}{G}_{1}^{0}-{\alpha }_{1,f}e\eta }{{kT}}\right)\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{g}_{-1}=\exp \left(-\displaystyle \frac{{\rm{\Delta }}{G}_{1}^{0}-{\rm{\Delta }}{G}_{{ad}}^{0}+{\alpha }_{1,b}e\eta }{{kT}}\right)\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{g}_{+2}=\exp \left(-\displaystyle \frac{{\rm{\Delta }}{G}_{2}^{0}-{\rm{\Delta }}{G}_{{ad}}^{0}-{\alpha }_{2,f}e\eta }{{kT}}\right)\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{g}_{-2}=\exp \left(-\displaystyle \frac{{\rm{\Delta }}{G}_{2}^{0}+(1-{\alpha }_{2,b})e\eta }{{kT}}\right)\end{eqnarray} \tag{ 12 }$$

where ${\rm{\Delta }}{G}_{i}^{0}$ are the activation energies of the assumed intermediate reactions, ${\rm{\Delta }}{G}_{{ad}}^{0}$ is the adsorption energy of the assumed intermediate species, α_i,f and α_i,b are the forward and backward transfer coefficients of reaction i respectively, and η is the reaction overpotential. It is worth noting that the kinetic model is likely a major oversimplification of the true reaction mechanism, with a number of more complex mechanisms proposed, ²¹ hence no insights regarding the nature of the OER are considered in this work. Furthermore, the original model assumed α_i,f + α_i,f = 1, however the use of fractional stoichometric coefficients in Eqs. 5 and 6 implies that the assumed reactions are not single electron transfer reactions. Therefore both forward and backward transfer coefficients are fitted separately which helped to more accurately reproduce the reaction kinetics.

The inclusion of the liquid saturation, s_L , in Eq. 7 accounts for blocking of the active area by gas bubbles and is the method used to induce mass transport losses. An alternative method is to use a concentration ratio and compute the concentration of water molecules at the reaction site, which allows for the possibility of water vapor within a gas bubble maintaining the OER at a reduced level. Reducing the active area was chosen to determine the worst case scenario with regard to the mass transport limitations predicted by the model. This method is commonly used in the literature however validation is required to determine the most appropriate method of accounting for mass transport limitations.

The hydrogen evolution reaction in the cathode proceeds on a platinum catalyst, which has been extensively researched in the literature. There are a number of models that can capture the kinetics of the reaction, ^65–67 however this work uses the model proposed by Wang et al. ²⁷ for the hydrogen oxidation reaction (HOR) which was applied by Elbert et al. ⁶⁸ to study both the HOR and HER on Pt/Ru core-shell catalysts. The model is given as:

$$\begin{eqnarray}&&j={A}_{v}{j}^{* }\left({g}_{+T}{\left(1-{\theta }_{H}\right)}^{2}-{g}_{-T}{\theta }_{H}^{2}+{g}_{+H}(1-{\theta }_{H})-{g}_{-H}{\theta }_{H}\right)\end{eqnarray} \tag{ 13 }$$

where the intermediate absorbed hydrogen coverage is denoted θ_H and the g_i terms are exponentials given as:

$$\begin{eqnarray}&&{g}_{+T}=\exp \left(-\displaystyle \frac{{\rm{\Delta }}{G}_{+T}^{* 0}}{{kT}}\right)\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{g}_{-T}=\exp \left(-\displaystyle \frac{{\rm{\Delta }}{G}_{+T}^{* 0}+{\rm{\Delta }}{G}_{{ad}}^{0}}{{kT}}\right)\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{g}_{+H}=\exp \left(-\displaystyle \frac{{\rm{\Delta }}{G}_{+H}^{* 0}-0.5\eta }{{kT}}\right)\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{g}_{-H}=\exp \left(-\displaystyle \frac{{\rm{\Delta }}{G}_{+H}^{* 0}+{\rm{\Delta }}{G}_{{ad}}^{0}+0.5\eta }{{kT}}\right)\end{eqnarray} \tag{ 17 }$$

where ${\rm{\Delta }}{G}_{i}^{* 0}$ are the activation energies of the intermediate reactions. The values used in this work are not fitted as the losses incurred from the HER are significantly smaller than other losses, such as the OER and ohmic losses. The activation energies reported by Wang et al. ²⁷ from rotating disk electrodes experiments are used in this work.

The pore size distribution model is used to characterize the porous media present in the catalyst and transport layers. The PSD is used to relate the pressures computed by the Darcy equations to the liquid saturation at each point in the domain, and allows for the permeability of the gas and liquid phases to be computed. The model conceptualises the porous medium as being a bundle-of-capillaries and is described in greater detail in references, ^54,55 with an updated formulation provided in reference ⁶⁰ which is used in this work. The gas and liquid pressure is used to compute the capillary pressure within each pore, which is used to compute the critical radius according to the Young-Laplace equation:

$$\begin{eqnarray}&&{r}_{c}=-\displaystyle \frac{2\gamma \cos \theta }{{p}_{c}}\end{eqnarray} \tag{ 18 }$$

where γ is the surface tension of the water, θ is the contact angle of the porous medium and the capillary pressure, p_c , is taken in this work as the difference between the liquid and gas pressures, i.e. p_c = p_l − p_g . The critical radius of a hydrophobic medium will be the smallest pore size that will be invaded by liquid water at a particular pressure, for a hydrophilic medium it will be the largest pore size invaded.

The PSD used in this work is given by

$$\begin{eqnarray}&&\displaystyle \frac{{dX}(r)}{{dr}}=\displaystyle \frac{{{dX}}_{{HI}}(r)}{{dr}}+\displaystyle \frac{{{dX}}_{{HO}}(r)}{{dr}}\end{eqnarray} \tag{ 19 }$$

where X(r) is the volume of the liquid invaded pore space, normalised to the total pore volume, i.e. X(r) = V_lw (r)/V_p . X(r) can be described using a log-normal distribution, which is given as:

$$\begin{eqnarray}&&\displaystyle \frac{{dX}(r)}{{dr}}=\displaystyle \sum _{i={HI},{HO}}{F}_{i}\displaystyle \sum _{k}\displaystyle \frac{{f}_{i,k}}{{{rs}}_{i,k}\sqrt{2\pi }}{E}_{i,k}^{(s)}(r)\end{eqnarray} \tag{ 20 }$$

where ${E}_{i,k}^{(s)}(r)$ is given by:

$$\begin{eqnarray}&&{E}_{i,k}^{(s)}(r)=\exp \left(-{\left[\displaystyle \frac{\mathrm{ln}r-\mathrm{ln}{r}_{i,k}}{{s}_{i,k}\sqrt{2}}\right]}^{2}\right)\end{eqnarray} \tag{ 21 }$$

The PSD can be comprised of a number of modes, k, each of which has a characteristic radius, r_i,k , and standard deviation, s_i,k , and take up a fraction, f_i,k , of the pore volume assigned to phase i. The hydrophobic/hydrophilic portions of the PSD are denoted by subscript i, with each having a fraction, F_i , of the total pore space, i.e., F_HO + F_HI = 1.

The liquid saturation is computed based on the local capillary pressure, as all the pores larger than the hydrophobic critical radius will be water filled, as will all pores smaller than the hydrophilic critical radius. As such the liquid saturation is computed by integrating over all the pores that have been invaded by water:

$$\begin{eqnarray}&&{s}_{L}={\int }_{0}^{{r}_{c},{HI}}\displaystyle \frac{{{dX}}_{{HI}}}{{dr}}{dr}+{\int }_{{r}_{c},{HO}}^{\infty }\displaystyle \frac{{{dX}}_{{HO}}}{{dr}}{dr}\end{eqnarray} \tag{ 22 }$$

The analytical solution to this integral is given as:

$$\begin{eqnarray}&&{s}_{L}=\displaystyle \sum _{i={HI},{HO}}{F}_{i}\displaystyle \sum _{k}\displaystyle \frac{{f}_{i,k}}{2}\left[1+{\xi }_{i}\mathrm{erf}\left(\displaystyle \frac{\mathrm{ln}({r}_{c,i})-\mathrm{ln}({r}_{i,k})}{{s}_{i,k}\sqrt{2}}\right)\right]\end{eqnarray} \tag{ 23 }$$

where ξ_HI = 1 and ξ_HO = − 1. The computed liquid saturation defines the volume of the pore space that is filled with water, i.e.:

$$\begin{eqnarray}&&{s}_{L}=\displaystyle \frac{V(\mathrm{water}\,\mathrm{filled}\,\mathrm{pores})}{V(\mathrm{total}\,\mathrm{pore}\,\mathrm{volume})}\end{eqnarray} \tag{ 24 }$$

while the gas saturation is given by s_G = 1 − s_L .

The permeability can be computed by comparing the velocity of a fluid passing through a cylindrical pore according to Hagen-Poiseuille equation to that computed by the Darcy equation: ⁶⁰

$$\begin{eqnarray}&&{k}_{{lw}}=\displaystyle \frac{{P}_{l-l}{\epsilon }_{p}}{16}\displaystyle \sum _{i={HI},{HO}}{F}_{i}\displaystyle \sum _{k}{f}_{i,k}{r}_{i,k}^{2}\exp (2{s}_{i,k}^{2}\left[1+{\xi }_{i,k}^{k}(r)\right])\end{eqnarray} \tag{ 25 }$$

where

$$\begin{eqnarray}&&{E}_{i,k}^{(k)}=\mathrm{erf}\left(\displaystyle \frac{\mathrm{ln}({r}_{c,i})-\mathrm{ln}({r}_{i,k})}{{s}_{i,k}\sqrt{2}}-{s}_{i,k}\sqrt{2}\right)\end{eqnarray} \tag{ 26 }$$

where P_l−l is the probability of a water-filled pore being in contact with another water-filled pore. The gas permeability can be found in a similar manner:

$$\begin{eqnarray}&&{k}_{g}=\displaystyle \frac{{P}_{g-g}{\epsilon }_{p}}{18}\displaystyle \sum _{i={HI},{HO}}{F}_{i}\displaystyle \sum _{k}{f}_{i,k}{r}_{i,k}^{2}\exp (2{s}_{i,k}^{2}\left[1-{\xi }_{i,k}^{k}(r)\right])\end{eqnarray} \tag{ 27 }$$

where P_g−g is the probability of two gas filled pores being in contact. The probabilities are based on determining the water/gas cross-sectional area per unit volume and are given as:

$$\begin{eqnarray}&&{P}_{l-l}={\lambda }_{{PSD}}{\left(\displaystyle \frac{{a}_{c}}{{a}_{{\max }}}\right)}^{2},{P}_{g-g}={\lambda }_{{PSD}}{\left(1-\displaystyle \frac{{a}_{c}}{{a}_{{\max }}}\right)}^{2}\end{eqnarray} \tag{ 28 }$$

where a_c is the water/gas cross-sectional area at a given capillary pressure, a_max is the total area for the entire pore volume and λ_PSD is the pore interconnectivity factor, which is used as a fitting factor.

The boundary conditions for the potential of the electronic phase are zero volts at the cathode GDL/land, and is the applied potential at the anode PTL/land, while no flux is applied elsewhere. A no flux boundary condition is applied for the protonic potential at all boundaries. Ohm's law for the protonic potential is not solved in the transport layers so the boundary condition is essentially applied at the CL/transport layer interface. The liquid pressure boundary condition is artificially set to a high negative pressure at the cathode channel to ensure that the liquid saturation is zero. The anode boundary condition is set to atmospheric pressure plus half the pressure drop down the single serpentine channel. The gas pressure is set to atmospheric in the cathode channel, while on the anode the gas pressure is set to that required to generate a bubble half the size of the channel according to the Young-Laplace equation, i.e.:

$$\begin{eqnarray}&&{p}_{g}={p}_{l}+\displaystyle \frac{2\gamma }{r}\end{eqnarray} \tag{ 29 }$$

The boundary conditions used in the anode electrode are greatly simplified, and could be improved to better reflect the two-phase regimes that are present in the anode channel under different operating current densities. The transition between bubbly, slug and annular flow has been well studied in the literature, however the use of the two-phase model in this work is primary concerned with the impact of the liquid saturation in the ACL and how the differing catalyst conductivities impact the two-phase flow within the ACL and the subsequent impact on the cell performance. A study was performed on the boundary condition used for the gas pressure in the anode, and while the liquid saturation in the PTL reduced significantly with increasing gas pressure at the boundary, the liquid saturation in the ACL was hardly affected. This is because of the weak coupling of two-phase effects between the PTL and ACL as is shown below.

The model can be used to compute the ohmic resistance of the cell, allowing for comparison with experimentally obtained resistances, such as the HFR. Using the ohmic heating method, ⁶⁹ the resistance to charge transport through a finite element can be computed using:

$$\begin{eqnarray}&&{R}_{{H}^{+}}^{{eff}}=\displaystyle \frac{1}{{i}^{2}A}\int {\sigma }_{m}{\rm{\nabla }}{\phi }_{m}\cdot {\rm{\nabla }}{\phi }_{m}{dV}\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&{R}_{{e}^{-}}^{{eff}}=\displaystyle \frac{1}{{i}^{2}A}\int {\sigma }_{s}{\rm{\nabla }}{\phi }_{s}\cdot {\rm{\nabla }}{\phi }_{s}{dV}\end{eqnarray} \tag{ 31 }$$

where i is the current density of the cell and A is the geometric area of the cell.

The model described was implemented in OpenFCST, an in-house, open-source software developed to model electrochemical systems. ⁷⁰ Second order Lagrange finite elements are used, with approximately 3000 cells, for a total of approximately 50,000 degrees of freedom. The governing equations are linearised using Newton's method, the resulting linear system is solved using UMFPACK. A tolerance of 10⁻⁸ was used for the residual of the Newton solver. A polarization curve up to a cell potential of 2 V took approximately ten minutes to obtain, using an Intel®Core^TMi7-8700K CPU at 3.7 GHz.

The PSD for the titanium sinter was fitted to mercury intrusion porosimetry data obtained in house, and its pore connectivity factor was adjusted to fit the permeability value reported by Jung et al. ⁷² The resulting relative gas and liquid permeabilities are shown in Fig. S2 in the Supporting Information (SI). The contact angle for the PTL was reported by Bromberger et al. ⁷⁴ as 56°. This value was not fitted as two phase transport in the PTL was not found to influence gas saturation in the ACL and therefore did not affect the cell performance, as discussed in the results section.

Obtaining the PSD of the Ir and IrOx ACLs used in this work would require performing mercury intrusion porosimetry (MIP) or focused ion beam scanning electron microscopy (FIBSEM) on the layers. The former is very challenging given the small thickness of the electrode, while the latter is expensive and time consuming. As the PSD was not fitted to experimental data presented in this work, the model is used to assess the impact of different PSD parameters on predicted performance of the cell. Several PSDs were used, representing extreme cases discussed in the literature. Three PSDs were available: a PSD with small pores used by Kosakian, ⁶⁰ a PSD from Hegge et al. ⁷⁵ with medium pores sizes, and a PSD from Wrubel et al. ⁵⁰ which had large pores. The base case uses the small PSD, while the other two PSDs are studied to provide a contrast. Each PSD will be applied in turn to both the Ir and IrOx based ACLs studied in this work, in order to understand the interplay between the reaction distribution and the ACL microstructure. The PSD from Hegge et al. ⁷⁵ was obtained by performing FIBSEM on an IrRuOx layer from a commercial CCM, with pore sizes on the order of a few hundreds of nanometres, as shown in Fig. 2a. Wrubel et al. ⁵⁰ used a similar numerical model to this work to study PTLs that had been chemically etched to obtained tunable features and used an intrusion curve based on a Pt/C layer. This intrusion curve was used in this work as the starting point to generate PSD under the assumption of an 80° contact angle, resulting in a PSD with a majority of pores close to a micron in size. The resulting PSD did not generate an intrusion curve that perfectly matched that from Wrubel et al., as their curve did not show zero liquid saturation. However, the resulting PSD was used as it was sufficiently different from the other two used in this work and allowed for the impact of having very large pore radii to be studied. It should be noted that a portion of this PSD has a characteristic pore radius of 6 μm, which is larger than the thickness of the layer. These pores are included as they are necessary to capture the intrusion curve used by Wrubel et al. despite perhaps being non-physical, though they could be considered to represent cracks through the ACL. The PSDs from Hegge et al. and Wrubel et al. have characteristic pore radii that do not appear to be present in surface SEM images of the catalyst layers used in this work; examples of these images can be seen in reference. ⁵⁸ As such, a PSD used by Kosakian ⁶⁰ based on FIBSEM analysis of a Pt/C CCL used in PEM fuel cells ⁷³ is also studied, as its characteristic pore radius is of the order of tens of nanometres.

Zoom In Zoom Out Reset image size

There is considerable uncertainty regarding the remaining parameters used in the PSD model, namely the contact angle of the hydrophilic and hydrophobic portions of the ACL, the fraction of the PSD that is hydrophilic, and the pores connectivity, i.e. the parameter λ_PSD in Eq. 28, which determines the permeability of the layer. As such, these parameters will be part of a parametric study to determine their impact on performance. As a base case, it is assumed that the hydrophobicity of the ACL is caused by the ionomer, in particular its PTFE backbone. A value of 107° is used based on contact angle measurements on PTFE films. ⁷⁶ Note that the assumption that the hydrophobicity of the ACL is caused by the PTFE backbone is not meant to imply the ionomer is overall hydrophobic, as Nafion thin films are thought to be hydrophilic in the presence of liquid water. ⁷⁷ In this work, only a small fraction of the ACL is assumed to be hydrophobic, with a contact angle that is chosen on the basis that the PTFE backbone is likely the only part of the ACL that is hydrophobic. The hydrophilic contact is assumed to be 80° for the base case as was used by Wu et al. ⁴⁹ in their numerical simulations. Wang et al. ⁴² measured a contact angle of 22°–53° depending on whether the ACL was conditioned or not, while contact angle measurements on iridium oxide films used in bioelectrodes suggest a contact angle closer to 10°. ^78,79 The values of the pore interconnectivity factor, λ_PSD , for the small and medium PSD is adjusted in order to obtain reasonable permeability values. Hegge et al. reported permeabilities of the order of 4 × 10⁻¹⁵ m², leading to a pore interconnectivity factor of 0.75 for the medium PSD. For the small PSD, the value of 4.1 reported by Kosakian ⁶⁰ is used. The hydrophilic portion for these two PSDs is assumed to be 90% on the basis of the parametric study that follows. Kosakian assumed that only the smallest pores were hydrophobic, hence the use of only one mode for the hydrophobic portion. It is unknown which pores are hydrophobic in the PSD from Hegge et al., so each mode was included in the hydrophobic portion. For the large PSD, the intrusion curve used by Wrubel et al. implied both a hydrophilic and hydrophobic portion, as the curve showed part of the layer is invaded for both negative and positive capillary pressures, leading to a hydrophilic fraction of 80%. Wrubel et al. did not report permeability values so the interconnectivity factor for their PSD was assumed to be the same as that used for Hegge et al. The pore interconnectivity was found to have no impact on the two phase behavior for this PSD. The assumption of a contact angle of 80° was made on the basis that smaller contact angles require larger pores in order to reproduce the intrusion curve used by Wrubel et al. For example, a contact angle of 40° results in an even large PSD, with the majority of pores having a characteristic pore radius of 3.7 microns. The relative gas and liquid permeabilities for each PSD using the parameters described are shown in Fig. S3 in the SI.

The membrane ionic conductivity is adjusted to match the membrane resistance reported by Mandal et al. ³⁰ for a liquid equilibrated N-211. Mandal et al. also reported the resistance of the cell hardware to be 19.7 mΩ.cm², which included the cathode GDL, which is modelled in this work, as well as two current collectors and bipolar plates, which are not included in the model. The GDL electronic conductivity was taken from Kosakian et al., ⁷¹ which results in a resistance of approximately 7.5 mΩ.cm² being computed using the ohmic heating method, see Eq. 30. The remaining 12.2 mΩ.cm² from the cell hardware is attributed to components not included in these simulations, so the experimental data was corrected to account for this resistance. This resistance was also removed from the experimentally obtained HFR. The PTL conductivity was adjusted to match the resistance reported by Mandal et al., ³⁰ who used a four point probe method to obtain the in-plane conductivity, which was assumed to match the through-plane conductivity due to isotropic nature of a sintered titanium PTL. The electronic and protonic conductivities of the cathode catalyst layer were again obtained from Kosakian et al. ⁷¹ These layers were not found to contribute significantly to the resistance of the cell.

The conductivities of the solid, σ_s , and ionomer, σ_m , phases in the ACL are dependent on whether an Ir or IrOx catalyst is used. The IrOx ACL electronic and protonic conductivities were based on the values reported by Mandal et al., ³⁰ using the data at a relative humidity (RH) of 100%. The Ir ACL conductivities were measured using the same methodology for this work. It should be noted that the hydrogen pump method used to measure the protonic conductivity of the Ir ACL may overpredict the layer conductivity. The Ir catalyst is active to the hydrogen reaction, which may result in the charge being transported through the electronic phase. ^30,80 However if this is the case, the trends described in this work will still be valid, as the electronic conductivity will still be higher than the protonic. As such the reaction distribution in the Ir layer will not be changed. The effective protonic conductivity at 100% RH, with an ionomer loading of 5%, was measured as 5.03 × 10⁻³ S cm⁻¹, while the effective electronic conductivity was 0.048 S cm⁻¹. These values, and those reported in reference ³⁰ for the IrOx, might not be representative of the ACL conductivities in the presence of liquid water, therefore they were slightly adjusted so that the HFR computed numerically would match that obtained experimentally (including the adjustment from the cell hardware resistance). The HFR from the numerical model is obtained by summing the electronic and protonic resistances of each layer computed from ohmic heating, besides the CCL and ACL. The contribution of an individual CL to the HFR, at low current densities, can be estimated as: ^28,29

$$\begin{eqnarray}&&{R}_{{HF}}=\displaystyle \frac{L}{{\sigma }_{m}^{{eff}}+{\sigma }_{s}^{{eff}}}\end{eqnarray} \tag{ 32 }$$

where L is the thickness of the layer, and σ^eff is either the effective electronic or protonic conductivity of the CL. Equation (32) shows that the computed resistance will be dominated by the more conductive phase, which was assumed to be the protonic phase for the IrOx ACL, and the electronic for the Ir ACL, based on the conductivity measurements from Mandal et al., ³⁰ therefore only the conductivity of the more dominant phase was adjusted to match the HFR. For the Ir ACL, the electronic conductivity was adjusted from the measured value of 0.048 S cm⁻¹ to 0.029 S cm⁻¹. It can be expected that the electronic conductivity of the ACL will be lower when liquid equilibrated compared to at 100% RH, as the swelling of the ionomer may cause loss of connectivity between catalyst particles. ⁸¹ The IrOx protonic conductivity was increased from a measured value of 0.011 S cm⁻¹ to 0.026 S cm⁻¹, which again is reasonable given that the ionomer will likely be more conductive when in contact with liquid water. Ideally the conductivities of these layers should be measured in a liquid equilibrated state to remove the uncertainty noted here.

The kinetic parameters were obtained by fitting to the experimental data at low current densities. Typically, the Tafel equation is used to describe the reaction kinetics of the OER and is obtained by fitting a straight line on a Tafel plot of the kinetic region of the experimentally obtained polarization curve, that has been corrected to account for the HFR. Using this method failed to provide a reasonable fit for either catalyst, as did fitting the kinetic model from Ma et al. ²⁴ instead of the Tafel equation, up to a current density of 0.5 A cm⁻². In order to obtain kinetic parameters that would result in a close fit to the experimental polarization curve, it was necessary to account for the resistance of the ACL while fitting the kinetic model. To do this, 2-D simulations were run with the ACL only, in which the kinetic parameters were fit to experimental polarization data that had the resistance of every other layer removed. The other layers were removed to minimize the computational requirements for the fitting procedure, which used a pattern search algorithm implemented in Dakota. ⁸² The fitting parameters were the free energies of activation, ${\rm{\Delta }}{G}_{i}^{0}$, free energy of adsorption ${\rm{\Delta }}{G}_{{ad}}^{0}$, and the forward and backward transfer cofficients for each intermediate step, α_i,f and α_i,b . The current densities fitted ranged from 3-500 mA cm⁻², smaller current density ranges did not produce as good a fit, which may be because the change in Tafel slope occurs at a high potential, circa 1.4-1.5 V, as is shown in Fig. S1 the SI. The resulting kinetic parameters are shown in Table V, and were used to obtain the numerical polarization curve shown in Fig. 3a. While the OER is a multi-step reaction which could exhibit a changing Tafel slope, and there is some experimental evidence for the change in Tafel slope in OER reaction kinetics, ^23,24 it is possible that the use of the kinetic model from Ma et al. is actually capturing a different effect such as the interfacial resistance between the ACL and PTL. The extremely close fit between the numerical and experimental data suggests that the additional voltage loss at least mimics a change in Tafel slope. The fact that the ACL resistance needed to be accounted for when fitting the kinetic parameters raises a number of questions about how kinetic parameters should be determined in PEMWE, and what constitutes the kinetic region. Finally, it should be noted that mass transport effects were not accounted for during the fitting of kinetic parameters, as relatively low current densities were used.

Zoom In Zoom Out Reset image size

Figure 3a shows the fit of the numerical model to the experimental data for both the Ir and IrOx catalysts. A very close match is obtained, with the parameters indicated in Tables II to VI. The reason for studying the two catalysts noted in this work is the difference in the effective conductivities of the two phases in the ACL. The Ir based ACL has an electronic conductivity that is almost ten times that of the protonic phase, whereas the electronic conductivity of the IrOx based ACL is three orders of magnitude lower than the protonic. An interesting consequence of the fact that the electronic conductivity is higher for one catalyst and lower for the other is that the OER will tend to occur at opposite sides of the ACL, as shown in Fig. 4. Overall, the numerical model is capable of reproducing the performance for both catalysts despite the very different behavior exhibited between them.

Zoom In Zoom Out Reset image size

The reaction distribution at 100 mA cm⁻² for the IrOx catalyst in Fig. 4 is extremely localised to one side, requiring a log scale to allow the distribution to be seen. A video of how the reaction distribution changes throughout the polarization curve without the log scale is available in the SI. Therefore, any oxygen bubbles created in the IrOx ACL will be produced close the PTL and will have an easy escape from the ACL. Figure 5 shows the gas, liquid and capillary pressures in the IrOx ACL at 4 A cm⁻², along with the liquid saturation. The gas pressure does not show a significant increase, particularly in comparison to the Ir ACL which is discussed later. The liquid pressure hardly changes at all. This suggests that the water consumed by the OER is rapidly replenished by the PTL. This was found to be the case even if the average PTL liquid saturation dropped to 5%. Using Faraday's law, the flux of liquid water required to maintain a current density of 4 A cm⁻² is of the order of 1 × 10⁻⁴ ml/(cm².sec), which the model suggests can be easily transported by the PTL since the liquid saturation in the ACL is quite uniform even at high current densities. Zhao et al. ¹⁷ showed that water can be transported rapidly through micro-scale cracks and at the interfaces between Ti particles of very hydrophilic PTLs, which may explain the PTLs ability to transport water even at very high gas saturations. As such the liquid saturation in the ACL does not change with increasing current density and there is little impact on the predicted polarization curves from mass transport losses, as shown in Fig. 3b.

Zoom In Zoom Out Reset image size

For the Ir catalyst, the bubbles will be produced close to the PEM, which will become more localised as the current density increases, requiring bubbles to traverse most of the ACL in order to escape to the PTL. That the gas bubbles must take a tortuous path to escape the Ir based ACL is highlighted in Fig. 6, which shows the gas, liquid and capillary pressures in the ACL at 4 A cm⁻², along with the liquid saturation. Very high gas pressures are generated at the ACL/PEM interface because of the small, hydrophilic pores that comprise the small PSD. It is unlikely that the ACL can withstand such high pressures and the ACL would likely be damaged by gas trying to escape the layer. Despite the high gas pressures in the ACL, the resulting ACL liquid saturation is not changed significantly. The pores are of the order of tens of nanometres which theoretically would require extremely high pressures to displace the water, despite the relatively high contact angle of 80°. Figure 2b shows the intrusion curves for each of the PSDs examined in this work. The small PSD requires capillary pressures of the order of 7 bar to achieve a liquid saturation of 0.5, which is likely not feasible. The liquid pressure is not significantly changed, as the hydrophilic nature of the ACL allows for water to be quickly replenished at the OER sites. As such, there is very little change in the cell performance when the liquid saturation dependent active area is used, as shown in Fig. 3b. It would appear, therefore, to be very difficult to induce significant mass transport limitations with the small PSD and the catalysts used in this work, which is also true for the other PSDs, as will be discussed later. At this point it is worth noting that the homogeneous nature of the model used in this work fails to capture the impact of discrete blockages. For example, it would not capture a large gas bubble being trapped at a particular location, e.g. in the PTL at the ACL/PTL interface, and blanking a large section of the ACL. That the pressures, and therefore liquid saturation, and ACL properties are averaged throughout the domain means that there is always a hydrophobic fraction of pores that can move gas bubbles out of the ACL. This study would therefore suggest that research into mass transport limitations should be concerned with determining the conditions under which these discrete blockages can occur. For example, if performing neutron or X-ray imaging of the gas buildup in a PTL during operation, the focus should be on reporting the size and frequency of gas bubbles close to the PTL/ACL interface, rather than the average gas saturation across the entire PTL.

Zoom In Zoom Out Reset image size

There is some evidence for the very high gas pressures in the Ir based ACL in the literature. Mandal et al. ⁵⁸ obtained polarization curves up to 4 A cm⁻² for both the Ir and IrOx catalysts used in this work, with a very similar cell construction. Post operation SEM imaging revealed the presence of a large number of small cracks throughout the surface of the Ir ACL that were not present prior to operation. The cracks had a volcanic appearance and may have been caused by the gas bubbles building significant pressure deep within the ACL before rupturing the layer to release this pressure. These volcano-like cracks did not appear on the IrOx based ACL, indicating that a pressure buildup did not occur during operation, as predicted by the numerical model in this work. It should be noted that the Ir based ACL used a Ti felt PTL on the anode, while the IrOx ACL used a sinter. Furthermore, Mandal et al. ⁵⁸ studied the impact of adding carbon particles, as a pore former, to the ACL ink before fabrication, with the intent of corroding these particles during operation, thereby altering the porosity and PSD of the layer. With the addition of the pore former to the Ir ACL, the volcano-like cracks did not appear after operation, indicating that the change in porosity and PSD induced by the pore former prevented the gas pressure from building up. The fact that the ACL could be ruptured by the gas pressure buildup suggests that an ACL with a small PSD will have its microstructure altered during operation, resulting in an increase in the mean pore size and permeability of the layer. As such, it may not longer be possible to build high enough gas pressures to displace the water from the pores, as shown will be shown in the case of the larger PSDs discussed below.

Other possible evidence for the build up of gas pressure in the ACL is from Trinke et al. ⁸³ who studied gas permeation through the membrane during PEMWE operation, using an Ir black catalyst in the ACL that likely has similar effective conductivities to the Umicore Ir black catalyst studied in this work. The authors noted that the rate of oxygen permeation through the PEM to the cathode increases with increasing current density. The authors attributed this increased flux to an oxygen supersaturation in the ACL; this supersaturation may have been caused by an increased oxygen pressure generated at the PEM/ACL interface, due to the nature of the OER distribution with an Ir black catalyst. Figure 7 shows that the maximum gas pressure generated in the ACL is dependent on the current density for the Ir catalyst, whereas the pressure changes very little for the IrOx catalyst. As such, it is possible that Trinke et al. would not have made their observation had they used a different catalyst.

Zoom In Zoom Out Reset image size

The PTL used in this work is a Ti sinter, with its PSD and intrusion curve shown in Fig. 2. Figure 8a shows the average gas saturation in the in-plane direction in front of the channel, and is plotted at a number of points in the through plane direction. An explanatory figure is provided in Fig. S4 in the SI that shows the locations of this averaging. There is a build up of gas within the PTL with increasing current density, in agreement with observations by De Angelis et al. ⁵ using X-ray tomography. It should be noted that Seweryn et al. ³⁴ and Zlobinski et al. ³⁵ reported that the PTL gas saturation was independent of the current density. De Angelis et al. suggested that the higher resolution of their equipment was needed to reveal the dependence. The rapid increase in the gas saturation moving from the channel into the PTL was not observed by De Angelis et al., and may be an artifact of the Dirichlet boundary conditions used at the channel for the gas pressure. Of particular interest, however, is that at the ACL/PTL interface, there is a significant jump in gas saturation, which was also noted by Garcia-Salaberri. ⁴⁸ This is not unexpected, as the gas and liquid pressures must be continuous. Therefore, there must be a jump in saturation between the two layers due to the significant differences between their PSDs and their intrusion curves. These differences may be indirectly observed in the works by Wang et al., ^41,43 who showed that bubbles emerging from the PTL tend to be larger and tends to occur more slowly than those emerging directly from the ACL. The smaller bubbles are likely a result of the smaller pores in the ACL, with the higher frequency of emergence from the ACL being impacted by the higher hydrophilicity. This result means that experimental observations of the gas saturation in the PTL may not be representative of the gas saturation in the ACL, and in fact the PTL gas saturation could be significantly higher than the ACL. As such, the importance of mass transport limitations may be overpredicted based on these experimental observations. It is also interesting to note that ACL gas saturation appears to be independent of the PTL gas saturation, as the ACL gas saturation is very similar at each current density, while the PTL gas saturation varies significantly. This is a result of the much smaller PSD in the ACL and would indicate that determining ACL PSD properties is vital for understanding the impact of mass transport limitations in PEMWE.

Zoom In Zoom Out Reset image size

Again it should be pointed out that the homogeneous nature of the model means that if the model predicts a gas saturation of 30% at the PTL/ACL interface at 4 A cm⁻², that also means that there is 70% liquid saturation at the same location which will provide enough water for the ACL to wick into the layer. It is more likely that the gas bubbles that comprise that 30% gas saturation are concentrated into a single capillary finger, in which a significant portion of the ACL is blanked from receiving water, an effect not captured by volume-averaged models such as the one presented. This process would likely induce more significant mass transport losses and would depend on the morphology of the PTL to determine the size of the area blanked by the gas, and on whether the liquid permeability of the ACL is high enough to allow water to be drawn to the blanketed area from other portions of the ACL in the in-plane direction. However, the ability of the ACL to transport water in the in-plane direction is likely determined by the PSD, hence further study on the structure of these layers is required. This idea implies that greater insights into these losses may be obtained from experimental imaging works that report localised gas saturations close to the PTL rather than averaged across the layer, e.g. using tomography. ⁵

The use of the PSD model allows for the ACL microstructure to be parametrised and its impact on performance to be assessed, which was achieved by studying the medium and large PSDs shown in Fig. 2a. These PSDs have significantly larger pore sizes than the small PSD, though it should be noted that there are also differences in the hydrophilic volume fraction and the pore interconnectivity factor. The impact of the liquid saturation dependent active area using the medium and large PSD is given in Fig. 9. The IrOx catalyst shows a very minor change in performance, though it is larger for the largest PSD. This is due to the hydrophilic volume fraction, which is 0.9 for the small and medium PSD, and 0.8 for the largest PSD. The Ir catalyst also shows very minor changes for the medium PSD, however a noticeable difference is seen for the large PSD, indicating that the model used in this work is capable of inducing mass transport losses for certain PSDs. Interestingly, mass transport effects can be observed throughout the polarization curve, rather than at high current densities only. The importance of the hydrophilic fraction can be shown by considering the intrusion curves. The maximum capillary pressures at 4 A cm⁻² induced by the Ir based ACLs are 280 kPa, 12 kPa and 1.3 kPa for the small, medium and large PSDs respectively, as seen for example in Fig. 6 for the case of the small PSD. Each of these capillary pressures are close to the "knee" for their respective intrusion curves, as shown in Fig. 2b. As such, the impact of the saturation based kinetics is controlled by the hydrophilic portion of the PSD and does not vary with current density.

Zoom In Zoom Out Reset image size

Figure 10 shows the average ACL liquid saturation and maximum gas pressures observed for both PSDs and catalysts. The medium PSD shows significantly lower pressures are induced in the ACL for the Ir, as compared to the small PSD in Fig. 7, because the larger PSDs are more permeable to gas. The liquid saturation for the IrOx catalyst remains at 0.9, while the Ir only has a small deviation at higher current densities. In general, the medium PSD results in liquid saturation and pressures profiles that are more similar to each other as compared to the small PSD when changing catalysts. This effect is clearly demonstrated for the large PSD, where pressure and liquid saturation profiles are identical for both catalysts, indicating that the mass transport behavior of the layer is no longer dictated by the reaction distribution. The very large pores present in the large PSD, as well as the reduced hydrophilic content, greatly improve the transport of gas out of the layer, alleviating the large gas pressures predicted with the small PSD. It is interesting to note that the volcano-like cracks observed by Mandal et al. ⁵⁸ in Ir based ACLs post operation are of the order of a few microns across, in a layer that appears to have very small pores. The volcano-like cracks, once formed, are of a similar size to the 6 μm pores present in the large PSD, and so likely significantly improve the permeability of the ACL, as implied by the results obtained for the large PSD. In addition, the large PSD shows a slightly stronger coupling between the PTL and the ACL. This is shown in Fig. 8b, where the gas saturation in the ACL for the large PSD shows some variation with current density and therefore PTL gas saturation. That the PTL exerts a greater influence can be expected, as the intrusion curves for the the large PSD is the closest to that of the Ti sinter compared to the other two PSDs, as shown in Fig. 2b.

Zoom In Zoom Out Reset image size

It is interesting to note that despite the ACL liquid saturation in the largest PSD being identical for both catalysts, the Ir ACL sees a larger drop in performance in Fig. 9b compared to the IrOx when mass transport effects are accounted for. At 4 A cm⁻², the Ir ACL has a mass transport induced voltage loss of 13.4 mV, compared to 6.9 mV for the IrOx. This discrepancy can be explained by examining the kinetics of the OER, and is demonstrated using the Tafel equation for simplicity. A multi-step model was used in this work, however the model will behave as a Tafel model at potentials that are sufficiently higher than the potential at which the change in Tafel slope occurs (approximately 1.4-1.5 V for parameters used in this work). The Tafel equation, including the liquid saturation dependent active area, is given as:

$$\begin{eqnarray}&&i={s}_{L}{A}_{v}{i}_{0}\exp \left(\displaystyle \frac{2.303\eta }{b}\right)\end{eqnarray} \tag{ 33 }$$

where b is the Tafel slope in Volts per decade, which can be extracted from Fig. S1 in the SI, and i₀ is the exchange current density. The equation can be rearranged to obtain the overpotential, η:

$$\begin{eqnarray}&&\eta =\displaystyle \frac{b}{2.303}\left(\mathrm{ln}(i)-\mathrm{ln}({s}_{L}{A}_{v}{i}_{0})\right)\end{eqnarray} \tag{ 34 }$$

The voltage loss induced by a reduction in the liquid saturation can be obtained by subtracting the overpotential obtained at a liquid saturation of one (i.e., no mass transport losses), from the overpotential obtained with a lower liquid saturation. The result is given as:

$$\begin{eqnarray}&&{\eta }_{{s}_{L}=1}-{\eta }_{{s}_{L}}=\displaystyle \frac{b}{2.303}\mathrm{ln}\left({s}_{L}\right)\end{eqnarray} \tag{ 35 }$$

Using the equation above, the voltage loss from mass transport can be computed for each catalyst. At 4 A cm⁻², both catalysts had an average liquid saturation of approximately 0.765. At high cell potentials, the Ir catalyst has a Tafel slope of 118.2 mV dec⁻¹, giving a voltage loss of 13.74 mV, while the IrOx layer has a Tafel slope of 59.25 mV dec⁻¹, giving a voltage loss of 6.89 mV. These values match very closely with those obtained from the difference in the polarization curve. This result highlights the importance of understanding the activity of the catalyst when studying mass transport losses. It should be noted however, that the analytical method described will not be as accurate if the liquid saturation distribution depends strongly on the reaction distribution, or there is a significant variation throughout the ACL.

Due to the uncertainty in obtaining certain PSD parameters, such as the hydrophilic portion, pore connectivity factor and the contact angle for the hydrophilic and hydrophobic portions, a parametric study was carried out on these parameters, using the three PSDs. Both catalysts were studied, however the IrOx did not show significant changes in behavior, so the results are not shown here. The study varied one parameter at a time, at a current density of 4 A cm⁻², where mass transport limitations should be most significant. The hydrophobic contact angle was not found to induce any change in behavior for either PSD or catalyst so is not shown. As the gas pressure will always be higher than the liquid pressure, the capillary pressure in the ACL will always be negative. As such, it is not possible for the water to invade the hydrophobic pores, so changing the contact angle of these pores cannot affect the mass transport behavior.

Figures 11a, 11d and 11g show that changing the hydrophilic contact angle can have a significant effect on the average liquid saturation and maximum gas pressure in the ACL, as well as the cell potential. For low contact angles, the liquid saturation is close to 0.9 or 0.8, which is the hydrophilic volume fraction of the ACL. At contact angles close to 90°, the liquid saturation drops dramatically for each PSD, as the pores are unable to hold onto the liquid water, and voltage losses increase. The largest PSD sees a very large decrease in ACL liquid saturation causing a large increase in the voltage required to reach 4 A cm⁻², as large pores with a high contact angle require very little pressure to invade. As such it may be possible to induce very significant mass transport losses with this homogeneous PSD model. However, it is not likely that the contact angle of the iridium based catalysts are close to 90°, particularly in light of the values reported by Cruz et al. ⁷⁸ and Chung et al. ⁷⁹ which were closer to 10°. Note that varying the IrOx hydrophilic contact angle also showed the same behavior for the liquid saturation as the Ir did, i.e. at very high contact angles, the liquid saturation dropped significantly. The gas pressure for the smallest PSD is strongly impacted by the contact angle, as lower angles prevent the gas from invading liquid filled pores, requiring high pressures to escape through hydrophobic section, which only makes up the remaining 10% of the volume. The gas pressure for the other two PSDs is less affected, due to their larger pore sizes which prevent gas pressure buildup.

Zoom In Zoom Out Reset image size

The change in the hydrophilic volume fraction induces an almost linear change in the liquid saturation, up to a volume fraction of approximately 0.9 for the small and medium PSDs, as the ACL liquid saturation essentially matches the hydrophilic volume fraction. As was discussed previously, the hydrophobic portion of the ACL is never invaded by the water, and provides an increasingly permeable pathway for the gas with decreasing hydrophilic volume fraction. This result suggests that an ACL with a hydrophobic portion of 50% means that 50% of the active area is unavailable for the OER, causing a significant reduction in the catalyst utilization. At very high HI volume fractions, the trend deviates from linearity for the smaller PSDs, as the reduced hydrophobic volume fraction reduces the permeability of the layer, requiring higher gas pressures to move the gas out. The increased gas pressure allows for more of the hydrophilic portion to be invaded by the gas, thereby reducing the liquid saturation. However, significantly reducing the liquid saturation again requires very high gas pressures.

Finally, the pore interconnectivity factor reduces the liquid saturation in the ACL, as it makes the ACL less permeable. However, there is a corresponding rise in the gas pressure, to values that would likely damage the ACL. The larger the PSD, the less impactful the interconnectivity factor becomes.