Modeling Overpotentials Related to Mass Transport Through Porous Transport Layers of PEM Water Electrolysis Cells

The adopted electrochemical model determines the overpotentials related to mass transport through PTLs, η_mt,PTL. Most commonly, mass transport overpotentials are calculated by taking the difference between the half-cell reduction potential E at the transport path's starting point and at its destination. The potential at the starting point corresponds to the hypothetical case in which transport resistances are zero. For the transport through PTLs, the potentials at their surfaces on the flow field and catalyst layer side are evaluated. Figure 1b depicts these surfaces (blue and orange) as well as associated definitions of model boundaries Γ. We sum the mass transport overpotentials of aPTL and cPTL to obtain the total overpotentials η_mt,PTL and calculate the required potentials via the Nernst equation:

$$\begin{eqnarray}\begin{array}{rcl}{\eta }_{\mathrm{mt},\mathrm{PTL}} & = & {\eta }_{\mathrm{mt},\mathrm{aPTL}}+{\eta }_{\mathrm{mt},\mathrm{cPTL}}\\ & = & \ {E}_{{ac}}-{E}_{{af}}+{E}_{{cc}}-{E}_{{cf}}\\ & = & \ \displaystyle \frac{R\,T}{4\,F}\left[\mathrm{ln}\left(\displaystyle \frac{{a}_{{ac}}^{{{\rm{O}}}_{2}}}{{\left({a}_{{ac}}^{{{\rm{H}}}_{2}{\rm{O}}}\right)}^{2}}\right)-\mathrm{ln}\left(\displaystyle \frac{{a}_{{af}}^{{{\rm{O}}}_{2}}}{{\left({a}_{{af}}^{{{\rm{H}}}_{2}{\rm{O}}}\right)}^{2}}\right)\right.\\ & & +\ \left.\mathrm{ln}\left({\left({a}_{{cc}}^{{{\rm{H}}}_{2}}\right)}^{2}\right)-\mathrm{ln}\left({\left({a}_{{cf}}^{{{\rm{H}}}_{2}}\right)}^{2}\right)\right]\\ & = & \ \displaystyle \frac{R\,T}{4\,F}\left[\mathrm{ln}\left(\displaystyle \frac{{a}_{{ac}}^{{{\rm{O}}}_{2}}}{{a}_{{af}}^{{{\rm{O}}}_{2}}}\right)+2\,\mathrm{ln}\left(\displaystyle \frac{{a}_{{cc}}^{{{\rm{H}}}_{2}}}{{a}_{{cf}}^{{{\rm{H}}}_{2}}}\right)-2\,\mathrm{ln}\left(\displaystyle \frac{{a}_{{ac}}^{{{\rm{H}}}_{2}{\rm{O}}}}{{a}_{{af}}^{{{\rm{H}}}_{2}{\rm{O}}}}\right)\right].\end{array}\end{eqnarray} \tag{ 1 }$$

The last row in Eq. 1 represents the most common grouping of the logarithmic arguments in the electrochemical community. The subscripts of the species' activities a indicate the model boundary Γ at which they are evaluated. Further quantities in the equation are the universal gas constant R, the temperature T and the Faraday constant F, which can also be found in the list of symbols at the end of this work. Protons' activities are not taken into account, as we assume that they depend only on the proton transport through the catalyst layers and membrane.

We assume that species can be present either in a liquid or gaseous phase and that phase compositions are following Henry's and Rault's laws. In this case, activities correspond to mole fractions,⁴²

$$\begin{eqnarray}&&{a}^{\kappa }={\chi }^{\kappa }=\displaystyle \frac{{n}^{\kappa }}{{n}_{\mathrm{total}}}=\displaystyle \frac{{c}^{\kappa }}{c},\kappa \in \{{{\rm{O}}}_{2},{{\rm{H}}}_{2},{{\rm{H}}}_{2}{\rm{O}}\},\end{eqnarray} \tag{ 2 }$$

where n denotes the amount (expressed in moles) of mixture constituents, i.e. species. We evaluate n over an arbitrary control volume V and obtain the shown notation in molar concentrations, in which c without a superscript denotes the concentration of all constituents regardless of their species. In the general two-phase case, the activities are calculated by

$$\begin{eqnarray}&&{a}^{\kappa }=\displaystyle \frac{{n}_{l}^{\kappa }+{n}_{g}^{\kappa }}{{n}_{\mathrm{total}}}=\displaystyle \frac{{S}_{l}\,{c}_{l}^{\kappa }+{S}_{g}\,{c}_{g}^{\kappa }}{{S}_{l}\,{c}_{l}+{S}_{g}\,{c}_{g}}.\end{eqnarray} \tag{ 3 }$$

The symbol S in above equation denotes the phase saturation and the subscripts l and g stand for the liquid and gaseous phase, respectively.

The concentrations in Eq. 3 are determined from the phase pressures p_g and p_l and the following assumptions. We require that both phases are at steady state conditions and neglect nucleation and related supersaturation effects,⁴³ such that thermodynamic equilibrium applies. It follows for the gas phase that the constituents' partial pressures are proportional to their molar fractions according to Rault's law,

$$\begin{eqnarray}&&{p}_{g}^{\kappa }={\chi }_{g}^{\kappa }\,{p}_{g}.\end{eqnarray} \tag{ 4 }$$

We derive the gas density ρ_g from the ideal gas law, neglecting real gas effects, as compressibility factors for oxygen and hydrogen are in the range of 0.978 to 1.001 at investigated operating conditions ${{\vartheta }}_{{\rm{o}}{\rm{p}}}\in [40{}^{\circ }{\rm{C}},60{}^{\circ }{\rm{C}}],\,{{p}}_{{\rm{o}}{\rm{p}}}\in [{{p}}_{{\rm{a}}{\rm{m}}{\rm{b}}{\rm{i}}{\rm{e}}{\rm{n}}{\rm{t}}},50\,{\rm{b}}{\rm{a}}{\rm{r}}]$.⁴⁴ If the water vapor's partial pressure is assumed equal to the equilibrium vapor pressure ${p}_{g}^{{{\rm{H}}}_{2}{\rm{O}}}={p}_{g,\mathrm{sat}}^{{{\rm{H}}}_{2}{\rm{O}}}$, which can be calculated with standard thermodynamic formulations,⁴⁵ all mole fractions ${\chi }_{g}^{\kappa }$ are known and the species' molar concentrations are obtained via

$$\begin{eqnarray}&&{c}_{g}^{\kappa }=\displaystyle \frac{{\chi }_{g}^{\kappa }\,{\rho }_{g}}{{M}_{g}},\end{eqnarray} \tag{ 5 }$$

where M_g is the gas phase's mean molar mass. The gas viscosity μ_g depends on the composition of the gas phase and is calculated according to the Wilke method.⁴⁴ The composition of the liquid phase is calculated using Henry's law,

$$\begin{eqnarray}&&{{\chi }}_{{l}}^{{\kappa }}=\displaystyle \frac{{{p}}_{{g}}^{{\kappa }}}{{{K}}_{{\rm{H}}}^{{\kappa }}},\,{\kappa }\in \{{{\rm{O}}}_{2},{{\rm{H}}}_{2}\}.\end{eqnarray} \tag{ 6 }$$

The Henry volatility K_H as well as the liquid viscosity μ_l and the liquid water concentration ${c}_{l}^{{{\rm{H}}}_{2}{\rm{O}}}$ are calculated according to standard thermodynamic formulations.^45,46 The dissolved species' influence on μ_l and ${c}_{l}^{{{\rm{H}}}_{2}{\rm{O}}}$ is neglected and the solute concentrations are obtained from

$$\begin{eqnarray}&&{{c}}_{{l}}^{{\kappa }}=\displaystyle \frac{{{\chi }}_{{l}}^{{\kappa }}}{{{\chi }}_{{l}}^{{{\rm{H}}}_{2}{\rm{O}}}}\,{{c}}_{{l}}^{{{\rm{H}}}_{2}{\rm{O}}},\,{\kappa }\in \{{{\rm{O}}}_{2},{{\rm{H}}}_{2}\}.\end{eqnarray} \tag{ 7 }$$

The saturations for Eq. 3 as well as the phase pressures are determined by the mass transport model. The model's domain definition and boundary conditions are chosen according to the experimental setup described by Suermann et al.¹³ The transport through anode and cathode PTLs is modeled within two distinct solution domains Σ_a and Σ_c, as seen in Fig. 1b. The flow field in the experimental setup is comprised of five 20 mm long, straight channels. In such a differential cell setup, gradients of any quantity are negligible along the channels. We therefore neglect transport inside PTLs in this longitudinal direction. Instead of modeling all five channels along the other in-plane direction, we assume the same conditions under each channel and model only one. We make use of symmetry axes by bordering the solution domain along them and defining no-flow conditions on them, as shown by the dashed lines in Fig. 1a. The oxygen and hydrogen evolution reactions (OER and HER, see definitions above) are modeled as mass flux boundary conditions at the interfaces between PTLs and catalyst layers (CL). This corresponds to the assumption that both reactions and nucleation processes take place at the boundaries Γ_ac and Γ_cc or inside the CLs. At the boundary Γ_ac, the following mass fluxes j^κ are prescribed according to Faraday's law:

$$\begin{eqnarray}&&{j}_{a}^{{{\rm{H}}}_{2}{\rm{O}}}=\left(1+2\,{n}_{d}\right)\displaystyle \frac{i\,{M}^{{{\rm{H}}}_{2}{\rm{O}}}}{2\,F}\quad \mathrm{and}\quad {j}^{{{\rm{O}}}_{2}}=-\displaystyle \frac{i\,{M}^{{{\rm{O}}}_{2}}}{4\,F}.\end{eqnarray} \tag{ 8 }$$

We denote with i the current density and the signs of fluxes correspond to the through-plane axis x (positive in the anode-to-cathode direction). In addition to the OER, water flux due to electroosmotic drag is also considered by introducing a non-zero drag coefficient,

$$\begin{eqnarray}&&{n}_{d}=0.0134\,\displaystyle \frac{T}{\,{\rm{K}}}+0.03.\end{eqnarray} \tag{ 9 }$$

Equation 9 describes the experimentally observed dependency of the drag coefficient on membrane temperature.⁴⁷ On the cathode side boundary Γ_cc, the following mass fluxes are defined:

$$\begin{eqnarray}&&{j}_{c}^{{{\rm{H}}}_{2}{\rm{O}}}=2\,{n}_{d}\displaystyle \frac{i\,{M}^{{{\rm{H}}}_{2}{\rm{O}}}}{2\,F}\,{\rm{and}}\,{j}^{{{\rm{H}}}_{2}}=\displaystyle \frac{i\,{M}^{{{\rm{H}}}_{2}}}{2\,F}.\end{eqnarray} \tag{ 10 }$$

Mass fluxes other than those in Eqs. 8 and 10, particularly those due to gas crossover, are neglected in the current model.

For Eqs. 8 and 10, an assumption for the current density distribution along the in-plane coordinate y is required. We are not aware of any published or unpublished current density mapping results that differentiate the states under channels and ribs. We have found, however, that varying current densities by ±12.5% around a mean value changes model results for averaged η_mt,PTL by less than ±1% compared to a case with constant current density i(y) = const. (results not shown). Due to the relatively small influence and the fact that no experimental results are available, we have assumed constant current densities for all shown simulations. It should be noted that in this work, we always evaluate η_mt,PTL overpotentials averaged over the y coordinate as we compare model results to integrally measured data without any spatial distribution.

On the domains' flow field sides, outflow boundary conditions are defined at the channels and no-flow boundary conditions are set under the ribs. We assume for the channels that the liquid pressure inside them is known and that all gas is immediately transported away from the PTL surface. The latter assumption is based on the aforementioned differential cell concept and the fact that convection was reported for both anode and cathode flow fields¹³ and the channel lengths of 20 mm are considered too short to facilitate significant gas accumulation. To this end, we define the following conditions on the boundaries Γ_af and Γ_cf: p_l = p_op and S_g = 0, as shown in Fig. 1b.

By defining a capillary pressure boundary condition p_c > 0 instead of the saturation condition (S_g = 0), it is possible to account for the case when the gas phase requires a higher pressure than p_l in order to exit the PTL. As the herein considered PTLs are hydrophilic, the capillary pressure is defined as p_c = p_g − p_l. The said phenomenon was observed in microfluidic experiments at so-called limiting throats.³² It is not known if limiting throats play a role for the herein investigated materials. If so, for the investigated PTL type with the product identification T10, the limiting throat diameter could be approximated as 25.9 μm. This diameter has been associated with bottleneck effects in mercury intrusion porosimetry data for T10.⁴⁸ A corresponding capillary pressure boundary condition can be calculated via the Young-Laplace equation,

$$\begin{eqnarray}&&{p}_{c}=4\,\displaystyle \frac{{\gamma }\,\cos ({\theta })}{d},\end{eqnarray} \tag{ 11 }$$

assuming an often reported contact angle of 70° for this kind of materials.^12,26,49 Applying this capillary pressure boundary condition on Γ_af and Γ_cf, however, leads to virtually the same model results as the saturation condition S_g = 0 (less than 1% difference in η_mt,PTL for i ≥ 0.1 A cm⁻² at default conditions ϑ_op = 50 °C and p_op = p_ambient, results not shown). S_g = 0 is therefore preferred for T10. The influence of limiting throats could not be estimated for the other investigated materials due to lack of data, which is further discussed in the results section of this work.

In addition to the domain definition and the boundary conditions, the model consists of the governing equation, which is a steady-state mass conservation of the following form:

$$\begin{eqnarray}&&\displaystyle \sum _{{\beta }}{\rm{\nabla }}\cdot \left({{j}}_{{\beta },{\rm{d}}{\rm{i}}{\rm{f}}}^{{\kappa }}+{{j}}_{{\beta },{\rm{a}}{\rm{d}}{\rm{v}}}^{{\kappa }}\right)=0,\,{\kappa }\in \{{{\rm{O}}}_{2},{{\rm{H}}}_{2},{{\rm{H}}}_{2}{\rm{O}}\},{\beta }\in \{{l},{g}\}.\end{eqnarray} \tag{ 12 }$$

The governing equation incorporates both component transport and fluid flow, which are assumed to follow Fick's and Darcy's laws:

$$\begin{eqnarray}&&{j}_{\beta ,\mathrm{dif}}^{\kappa }=-{M}^{\kappa }\,{D}_{\mathrm{pm},\beta }\,{\rm{\nabla }}{c}_{\beta }^{\kappa },\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{j}_{\beta ,\mathrm{adv}}^{\kappa }=-{M}^{\kappa }\,{c}_{\beta }^{\kappa }\,\displaystyle \frac{{k}_{r\beta }}{{\mu }_{\beta }}\,K\,{\rm{\nabla }}{p}_{\beta }.\end{eqnarray} \tag{ 14 }$$

In Eq. 13, the effective diffusion coefficient D_pm,β is calculated according to the parallel pore model⁵⁰ and a Millington-Quirk-type tortuosity τ⁵¹:

$$\begin{eqnarray}&&{{D}}_{{\rm{p}}{\rm{m}},{\beta }}={\phi }\,{{\tau }}_{{\beta }}\,{{S}}_{{\beta }}\,{{D}}_{{\beta }},{{\tau }}_{{\beta }}={{\phi }}^{1/3}\,{{S}}_{{\beta }}^{7/3},\,{\beta }\in \{{l},{g}\}.\end{eqnarray} \tag{ 15 }$$

In this equation, ϕ is the PTL's porosity and the phases' diffusion coefficients D are calculated according to Ref. 44. For the material T10, diffusion of dissolved oxygen and hydrogen is negligible at current densities i ≥ 1 mA cm⁻² (diffusion's contribution to ${j}^{{{\rm{O}}}_{2}}$ and ${j}^{{{\rm{H}}}_{2}}$ is less than 0.05% at default conditions ϑ_op = 50 °C, p_op = p_ambient, as shown in the appendix), whereas vapor diffusion might play a role at high gas saturations, cf. Eq. 15.

In Eq. 14, the intrinsic and relative permeabilities are K and k_r, respectively. The underlying continuum approach behind Eq. 14 assumes each phase to be continuously distributed throughout the whole solution domain.⁴¹ It follows from this assumption that the two-phase flow can be parametrized with saturation dependent functions for the capillary pressure p_c(S_l) and for the relative permeabilities k_rl(S_l) and k_rg(S_l). These parameter functions are assumed to be unique and monotonic. For this reason, the model yields monotonic spatial profiles for both phase saturations and phase pressures, as seen in Fig. 1c, with pressure gradients $\parallel {\rm{\nabla }}{p}_{\beta }\parallel$ being generally higher at locations with low phase saturation S_β. Note that in this work, we focus on cases with balanced operating pressures such that p_l is identical on Γ_af and Γ_cf. Accordingly, there is a pressure jump between both sides of the membrane that corresponds to the electroosmotic pressure differential, as shown in Fig. 1c. In this work, the pressure differential results from the PTL models and is not calculated by an electroosmosis model.

Including potential energy in Eq. 14 did not influence simulation results by more than 0.6% at default conditions (ϑ_op = 50 °C, p_op = p_ambient, results not shown), therefore gravitational effects are neglected in this work. Furthermore, model results are calculated under the assumption of isothermal conditions. This simplification has limited effect on the model results, as discussed in the results section.

To summarize, the presented mass transport model takes the current density i, the operating pressure p_op and the operating temperature ϑ_op as inputs and also needs parameters for the intrinsic permeability K, the porosity ϕ, the capillary pressure function p_c(S_l) and the relative permeabilities k_rβ(S_l). The model's outputs are the phase saturations and phase pressures, the latter of which determine the species concentrations ${c}_{\beta }^{\kappa }$.

The mass transport model is solved using the DuMuX framework.⁵² The solution of the governing Eq. 12 is found by time integration of its transient form,

$$\begin{eqnarray}&&\displaystyle \sum _{\beta }\phi \,{M}^{\kappa }\,\displaystyle \frac{\partial ({S}_{\beta }\,{c}_{\beta }^{\kappa })}{\partial \,t}+\displaystyle \sum _{\beta }{\rm{\nabla }}\cdot \left({j}_{\beta ,\mathrm{dif}}^{\kappa }+{j}_{\beta ,\mathrm{adv}}^{\kappa }\right)=0,\end{eqnarray} \tag{ 16 }$$

from a fully liquid saturated initial state (S_l(t₀) = 1 over Σ_a and Σ_c) until steady-state conditions are reached. The box method is adopted for spatial discretization and a backward Euler scheme is used for time integration.⁵⁰ A grid resolution of 80 by 160 elements was employed after performing a grid sensitivity study (based on <1% change in modeled overpotentials when doubling grid resolution in both dimensions). The solution process, furthermore, includes Newton-Raphson iterations in order to treat non-linearities and the coupled nature of solved equations.

The examined PTL materials are three sintered titanium powders from type SIKA T5, T10 and T20 produced by GKN (DE). As the materials were initially manufactured for the filter industry, the numbering corresponds to the maximum particle diameter (in  μm) passing the material in through-plane direction. All materials have a nominal thickness of 1 mm and have been analyzed with X-ray tomographic microscopy by Suermann et al.¹³ The porosities, intrinsic permeabilities and pore size distributions have been derived by flow computations and geometric processing using tomography data. The intrinsic permeabilities have been found to be isotropic.¹³ The main material properties are listed in Table I.

Table I.  Porous media properties of investigated PTL materials according to Ref. 13.

	T5	T10	T20
Permeability K/(10−12 m2)	1.4 ± 0.0	3.4 ± 0.1	7.7 ± 0.2
Thermal conduct. σ/(W m−1 K−1)	11.6 ± 0.3	9.6 ± 0.4	11.0 ± 0.7
Porosity ϕ/10−2	30 ± 2	35 ± 2	33 ± 2
Mean pore $d/\mu {\rm{m}}$	30	38	57

For T10's capillary pressure function p_c(S_l), Lettenmeier et al.⁶ have calculated curve points with the pore morphology method (material named "sintered Ti" therein). We adopt this curve and interpolate between points using cubic hermite splines, see Fig. 2a. For the materials T5 and T20, no such calculations are available and their capillary pressure curves can only be estimated, as discussed later in this work.

Zoom In Zoom Out Reset image size

Pore scale calculations of the relative permeabilities are not available for the investigated materials. We therefore apply the measured relative gas permeability curve of a sintered powder PTL, which is similar to the herein investigated materials (see material "sample 3" in Ref. 33). The k_rg values, referred to in the publication as "normalized dry cross-sectional area," were obtained by capillary flow porosimetry with a quasi-static liquid phase. Global saturations were determined by weighing the sample. Global liquid saturations were measured in the range from 0 to 0.75 (upper limit is determined by breakthrough) which are herein linearly scaled to the [0, 1] range in order to obtain local saturations (local liquid saturations at throats, which dictate k_rg, are generally higher than global saturations, cf.³²).

We are not aware of any published or unpublished measurements of liquid permeabilities for sintered powder PTLs and apply a curve obtained with the Mualem model.⁵³ The model is standard in soil sciences and constructs the liquid permeability based on the capillary pressure curve:

$$\begin{eqnarray}&&{k}_{{rl}}\left({S}_{l}\right)=\sqrt{{S}_{l}}{\left({\int }_{0}^{{S}_{l}}\displaystyle \frac{{\rm{d}}{S}_{l}}{\left.{p}_{c}({S}_{l}\right)}\right)}^{2}{\left({\int }_{0}^{1}\displaystyle \frac{{\rm{d}}{S}_{l}}{\left.{p}_{c}({S}_{l}\right)}\right)}^{-2}.\end{eqnarray} \tag{ 17 }$$

Both an alternative standard model by Burdine⁵⁴ and a constant liquid permeability k_rl(S_l) = 1 yield virtually the same results as the Mualem model (difference in calculated overpotentials is less than 0.01% at default conditions ϑ_op = 50 °C, p_op = p_ambient, i ≤ 4 A cm⁻², results not shown) and consequently, the model result are found to be negligibly sensitive to k_rl.

The considered PTL materials have also been characterized by operando measurements and a cell voltage breakdown.¹³ The voltage breakdown yields the overpotential η_mtx which includes all types of mass transport processes, cf. Fig. 3a.

$$\begin{eqnarray}&&{\eta }_{\mathrm{mtx}}={E}_{\mathrm{cell}}-{E}_{\mathrm{th}}^{0}-{\eta }_{\mathrm{act}}-{\eta }_{{\rm{\Omega }}}.\end{eqnarray} \tag{ 18 }$$

Zoom In Zoom Out Reset image size

Here, E_cell represents the measured cell voltage, ${E}_{\mathrm{th}}^{0}$ the thermodynamic cell voltage, η_act the activation overpotential and η_Ω the ohmic overpotential. The ohmic overpotentials are obtained from high frequency resistance (HFR) measurements and activation overpotentials are determined via the Butler-Volmer equation considering the transition state theory and the concept of a rate determining step.⁵⁶ The η_mtx overpotentials include a number of different loss types, where anodic and cathodic contributions are merged. Such types of losses contain, but are not limited to, mass transport resistances to and from the electrochemical active surface, proton transport resistances within the CL and current density dependent hydrogen and oxygen gas crossover through the membrane. The gas crossover leads to a concentration increase in the ionomer of the CL shifting the thermodynamic equilibrium cell voltage.⁵⁷ The CL's proton transport resistances are calculated using further experimental data with identical materials using a transmission line model approach.⁵⁵ According to the presented Nernst-type electrochemical model, proton transport is isolated from considered mass transport overpotentials:

$$\begin{eqnarray}&&{\eta }_{\mathrm{mt}}={\eta }_{\mathrm{mtx}}-{\eta }_{{{\rm{H}}}^{+},\mathrm{CL}}.\end{eqnarray} \tag{ 19 }$$

See Table II for an overview of overpotentials.

Table II.  Overview of examined overpotential types; in this work, experimental data is available for η_mtx, ${\eta }_{{{\rm{H}}}^{+},\mathrm{CL}}$, η_mt and η_mt,p−dep (as a result of voltage breakdown, mass transport overpotentials may also contain overpotentials due to other unknown processes).

Overpotential	Related mass transport process
ηmtx	All
${\eta }_{{{\rm{H}}}^{+},\mathrm{CL}}$	CL proton transport
ηmt	All excluding CL proton transport
ηmt,p−ind	pop-independent part of ηmt
ηmt,p−dep	pop-dependent share of ηmt
ηmt,PTL	Exclusively through bulk PTL

Variations of operating pressures showed a decrease of mass transport overpotentials with increasing operating pressure, especially for pressure configurations where both the anode and the oxygen sides are pressurized.⁵⁸ Moreover, η_mt overpotentials converge in the higher two-digit pressure range (measured in bar at symmetric pressure conditions) for all investigated PTL types,¹³ as shown in Fig. 3b. This convergence allows a division of η_mt overpotentials into pressure-dependent (above the converged value) and pressure-independent (equal to the converged value) overpotentials: η_mt,p−dep and η_mt,p−ind, respectively. The convergence value is found to be approximately the same for all investigated materials. Thus, a PTL-independent convergence value is found. From this, we derive the hypothesis that all PTL-related mass transport overpotentials η_mt,PTL are part of pressure-dependent mass transport overpotentials η_mt,p−dep. In other words, we assume that mass transport through PTLs does not affect η_mt,p−ind overpotentials, which account, depending on current density and PTL material, for ca. 40% to 75% of η_mt overpotentials at intermediate and high current densities (i.e. at 0.5 A cm⁻² ≤ i ≤ 4 A cm⁻² and default operating conditions ϑ_op = 50 °C, p_op = p_ambient).¹³ Model results for PTL-related overpotentials are therefore compared to the experimental data on η_mt,p−dep in this work. Breakdown data on η_mt,p−dep is available for different operating conditions (current density, pressure) as well as for different PTL types (T5, T10, T20). Although this already constitutes a comparatively broad data set, supplementary measurements would be helpful for further breakdown of η_mtx overpotentials, e.g. at different temperatures and pressures, differential pressure conditions (in both directions) or using a reference electrode separating anodic and cathodic contributions.

The adopted continuum model with previously described parameters yields the results for anode- and cathode-specific bulk PTL mass transport overpotentials, η_mt,aPTL and η_mt,cPTL, that are shown in Fig. 4. On both anode and cathode side, the PTL consists of the T10 material and operating conditions correspond to ambient pressure and an operating temperature of 50 °C. The overpotentials attributed to the cathode PTL are found to contribute by 68% to 71% to the total PTL-related mass transport overpotentials η_mt,PTL.

Zoom In Zoom Out Reset image size

Both model results and experimental results show that overpotentials rise with current density due to increasing mass fluxes, as shown in Fig. 4. As previously described, oxygen and hydrogen are transported mainly by advection, as diffusive transport is negligible at considered conditions. The increasing gas flow results in higher gas saturations and gas pressures on the PTLs' catalyst layer side (opposite to the flow field). This is also evident from saturation and pressure profiles along the through-plane coordinate x, Fig. 5. The pressurization of gas additionally leads to heightened oxygen and hydrogen concentrations in both liquid and gas phase, cf. Eqs. 4 and 6. All of these effects contribute to a gain in oxygen and hydrogen activities, a decrease in water activities and consequently to an increase in mass transport overpotentials. It has been observed for modeled overpotentials that the changes in saturations are more decisive than relative changes in concentrations, cf. Fig. 1 in supplementary data (is available online at stacks.iop.org/JES/167/114511/mmedia). Note that at considered current densities, gas pressure drops are two to three orders of magnitude higher than liquid pressure drops, as seen in Fig. 5. This is due to a higher resistance for gas flow, which results from low gas permeabilities k_rg at high liquid saturations, as shown in Fig. 2b. Gas pressure profiles are similar on the anode and the cathode side due to three facts: (1) molar concentrations ${c}_{g}^{{{\rm{O}}}_{2}}$ and ${c}_{g}^{{{\rm{H}}}_{2}}$ are the same at a given gas pressure (follows directly from the ideal gas law and identical equilibrium vapor pressures), (2) hydrogen's stoichiometric factor in water electrolysis is two times higher than oxygens's, (3) vapor saturated oxygen has a 1.93 times higher dynamic viscosity than vapor saturated hydrogen at gas pressures close to p_ambient.⁴⁴ It follows from Eqs. 8, 10 and 14 that the second and third fact nearly balance each other out for gas flow.

Zoom In Zoom Out Reset image size

The model and experimental results are in excellent qualitative and reasonable quantitative agreement, especially for current densities above 0.5 A cm⁻². As a quick reminder, the model results shown in Fig. 4 are generated without any calibration of parameters. However, as expected, the agreement is not as good at lower current densities, $i\lt 0.5$ A cm⁻², as briefly explained below. While the model considers transport losses as soon as Faradaic reactions start, i.e. $i\gt 0$ A cm⁻², the electrochemical breakdown method by definition ignores any transport effects below approximately 0.1 A cm⁻². An in-depth discussion of these different assumptions is presented later in this work.

So far we have mainly discussed the influence of the operating parameter i, which is the current density. From experimental findings we know that the mass transport overpotentials decrease with increasing operating pressure. From a theoretical, more fundamental perspective, fluid states change at high operating pressures mainly due to the increase in phase viscosities, dissolution and gas density, the first of which raises transport resistance while the latter two reduce gas saturations considerably. As gas saturations have a high influence on mass transport overpotentials, the model matches the experimentally observed falling trend in overpotentials with rising operating pressures, as shown in Fig. 6. Between model results and breakdown data, there is a difference considering the rate with which overpotentials decrease at rising p_op (as can be seen for the 10 bar case), yet a quantitative agreement is observed again at 50 bar balanced operating pressure.

Zoom In Zoom Out Reset image size

Another crucial operating parameter for the PEMWE is the temperature, as it influences all main overpotentials, as well as the fluid states and fluid properties. Model results at increasing operating temperatures are determined by gas expansion, increase in gas viscosity and outgassing of solutes. As a result of these effects, higher gas saturations and mass transport overpotentials are obtained at higher temperatures, as shown in Fig. 7. There is unfortunately no breakdown data available for comparison.

Zoom In Zoom Out Reset image size

Nevertheless, it is known that temperature gradients are present in PEMWE cells and temperature variations of about 8 K have been measured by Suermann et al.¹³ at i = 4 A cm⁻² and nominal operating temperature of 50 °C using the same PTL material T10 (with a thermal conductivity of about 9.6 W m⁻¹ K⁻¹, cf. Table I). The used model neglects temperature variations and assumes the nominal operating temperature everywhere in the solution domain. An incorporation of temperature gradients in the model would yield larger overpotentials compared to the isothermal results shown in Fig. 7. However, non-isothermal solutions for ϑ_op = 50 °C are expected to be still lower than the (brown) isothermal 60 °C curve shown in Fig. 7 and the model error due to the isothermal assumption is therefore less than 5 mV for current densities under 4 A cm⁻².

Briefly summarized, the results in this section illustrate that a very good match between modeled overpotentials and experimental data can be reached at different current densities and balanced operating pressures. In the presented case, where material properties were either known or could be approximated based on measurements, no parameter fitting was needed. Regarding PTLs' influence on mass transport overpotentials, both model and experimental results suggest that the transport through PTLs contributes to a pressure-dependent share of the mass transport overpotentials, which converge towards negligible values at pressures over ca. 80 bar. Furthermore, cathode-side mass transport through PTL is found to cause more than twice as much overpotentials as the anode-side mass transport processes if the same PTL material is used on the anode and the cathode side. The model calculations predict that PTL-specific mass transport overpotentials rise with increasing temperature mainly due to gas expansion, although experimental validation is still missing.

The PTL's fluid mechanical properties affect fluid flow and consequently advective transport of mass, which is found determinant for PTL-related mass transport overpotentials. Accordingly, the adopted model is sensitive to all material parameters considered in Darcy's law: (1) the intrinsic permeability K, (2) the relative permeabilities k_rl and k_rg and (3) the capillary pressure curve p_c(S_l). These parameters can be obtained either by an experimental approach, e.g. capillary flow porosimetry and saturation measurements,³³ or by conducting pore scale or pore network simulations.²⁸ If the saturation-dependent parameters (2) and (3) are not or only partially available, parameter models can be used to approximate these parameters. This is carried out below for the PTL material T10, which serves here as a reference.

The modified Kovács model⁵⁹ can be applied for predictions of the capillary pressure curve. This parameter model requires a known contact angle θ, which is not available for T10. However, the modified Kovács model yields a similar curve as the pore morphology calculation if θ = 70° is assumed according to literature data for similar materials,^12,26,49 as shown in Fig. 8a. Correspondingly, model results for η_mt,PTL are obtained with the modified Kovács model that agree qualitatively with experimentally determined overpotentials, yet overestimate them by up to 12 mV, as seen in Fig. 9.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The Mualem model⁵³ from soil sciences has been previously introduced for predicting relative permeabilities. Regarding k_rl(S_l), the model results were found insensitive at considered current densities i ≤ 4 A cm⁻². The predicted curve for k_rg(S_l), exhibits relatively small values at high liquid saturations S_l > 0.9, as shown in Fig. 8b. Hence, applying a Mualem-type gas permeability for the model yields large η_mt,PTL overpotentials at small current densities (Fig. 9), which overestimate breakdown data by ca. 20 mV. Using an alternative standard model for relative permeabilities, the Burdine model,⁵⁴ results in even (up to 17 mV) larger overpotentials compared to the Mualem model (results not shown). These findings illustrate that model results are highly sensitive to the choice of relative gas permeability, especially at low gas saturations, at which standard parameter models possibly underestimate gas permeabilities.

As mentioned in the introduction, the existence of an optimum pore size of the PTL has been reported in various experimental works regarding the reduction of total overpotentials. For a model-based optimization of the pore size, it is essential to understand how pore diameters affects the parameters to which the model is sensitive. As aforementioned, the adopted continuum model's results mainly depend on the parameters intrinsic permeability K, capillary pressure curve p_c(S_l) and relative gas permeability k_rg(S_l) at investigated current densities i ≤ 4 A cm⁻². In the following, the dependence of K and p_c(S_l) on the mean pore size and how model results are affected by changes in these parameters are discussed. The discussion of the relative gas permeability k_rg(S_l) is omitted here, as standard models, which could describe the pore size's influence, were found to insufficiently approximate k_rg(S_l) at large liquid saturations (at S_l > 0.9, see previous section).

We consider variations of the mean pore size that correspond to the previously introduced T5 and T20 materials: a 20% reduction and a 50% increase, respectively (${d}_{-20 \% }={d}_{{\rm{T}}5}=30\,\mu {\rm{m}}$ and ${d}_{+50 \% }={d}_{{\rm{T}}20}=57\,\mu {\rm{m}}$, cf. Table I). In a previous study,¹³ pore scale simulations yielded intrinsic permeabilites for T5 and T20 that follow a rising trend with increasing pore size, as also shown in Table I. This trend is reflected by the often-used Kozeny-Carman approximtaion,⁶⁰

$$\begin{eqnarray}&&K\propto {r}_{h}^{2},\end{eqnarray} \tag{ 20 }$$

in which the mean hydraulic radius r_h scales with the mean pore size d. Figure 10 shows T10's capillary pressure curve that is transferred to the varied pore sizes by applying the factor ${d}_{{\rm{T}}10}/d$ according to the Young-Laplace Eq. 11. It is evident from the above descriptions that reducing pore sizes corresponds with lowering the intrinsic permeability and increasing capillary pressures. While low intrinsic permeabilities lead to increased transport resistances, higher values for p_c lead to lower mass transport overpotentials due to the gas pressurization, as discussed above in context of the operating pressure. Increasing p_c values further facilitates gas flow, as slopes in the p_c(S_l) curve usually also become steeper, cf. Fig. 10, yielding higher gas pressure gradients $\parallel {\rm{\nabla }}{p}_{g}\parallel$ for a given saturation field. For the investigated cases, it is found that the influence of intrinsic permeability on the mass transports outweighs that of capillary pressure, as seen in Fig. 11. In general terms, however, the opposing effects of K(d) and p_c(d) support the experimental observations that an optimum pore size d exists.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A reduction of T10's mean pore size by 20% yields higher modeled η_mt,PTL values compared to T10's overpotentials, as shown in Fig. 11. This is in qualitative agreement with the corresponding breakdown results of T5 materials. Even though the model results for reduced pore size were obtained with T5's mean pore size and intrinsic permeability, a quantitative discrepancy of up to 30 mV can be observed regarding the experimental data (cf. blue curves in Fig. 11). The model results with T20's mean pore size and intrinsic permeability do not qualitatively reflect the experimental data, as T20's modeled overpotentials are smaller than those for T10 (due to the relatively high intrinsic permeability), whereas T20's experimental curve exhibits the largest of all overpotentials. The following three reasons can be given for the obvious discrepancy, especially for T20. First, assumed parameters might deviate from T5's and T20's actual ones, e.g. regarding the shape of saturation dependent parameter curves. The transfer of parameters between materials could also be complicated by the existence of discrete transport effects, such as limiting throats. It should be noted that all model results were obtained assuming the same contact angle for all materials (otherwise T10's capillary curve would have been transferred by the factor $\left({d}_{{\rm{T}}10}\,\cos (\theta )\right)/\left(d\,\cos ({\theta }_{{\rm{T}}10})\right)$ for T5 and T20, cf Eq. 11). Second, the adopted model assumes that reaction sites are uniformly distributed and very close to the PTLs' surface. As T20 consists of almost twice as large titanium particles as T5 and T10,¹³ it exhibits a comparably rough surface topology and its electrochemically active surface area deviates from the others.⁵⁵ Transport processes close to T20's (catalyst layer-side) interface might contribute to the high mass transport overpotentials that were experimentally determined for T20. Third, the experimental data underlies inaccuracies that originate in assumptions made in the breakdown method. These assumptions are discussed separately in the next section.

Regarding modeled η_mt,PTL overpotentials for the T10 material, the cathode PTL's share is in the range of 68% to 71%, as seen in Fig. 4. This is qualitatively and quantitatively in good agreement with electrochemical measurements using differential pressure operation (p_op,a < p_op,c, instead of balanced pressure operation p_op,a = p_op,c) with up to 30 bar⁶¹ and up to 50 bar⁶² differential pressure. However, further measurements at even higher pressures and in combination with a third (reference) electrode could yield a more precise data set concerning anodic and cathodic shares.

Nevertheless, all shown simulations consistently overestimate breakdown data at current densities below ca. 0.5 A cm⁻² due to a (fundamental) difference in model and breakdown assumptions. The applied breakdown method attributes all iR-free overpotentials (corrected by the ohmic overpotential) under ca. 0.1 A cm⁻² to activation overpotentials and thus mass transport overpotentials are neglected by definition. This breakdown method considers the Butler-Volmer approach, transition state theory and concept of a rate determining step⁵⁶ and ultimately results in the Tafel approach. Here, a Tafel-fit is made at small current densities in the linear region (typically 0.01 A cm⁻² ≤ i ≤ 0.1 A cm⁻²) of a Tafel-plot (iR-free E_cell versus log(i)) and extrapolated towards higher current densities. For more details of this overpotential breakdown method, we refer to the Ref. 63. The breakdown method's assumption is contrasted by the results of the adopted mass transport model, which show that oxygen and hydrogen are transported predominantly in the gas phase with considerable transport resistance at i = 0.1 A cm⁻². This model result is supported by operando gas saturation measurements in various anode PTLs (including T5) yielding 8% to 46% gas saturations at i = 0.1 A cm⁻².^39,40 Consequently, bubble nucleation already occurs at those relatively small current densities causing mass transport resistances. The current density at which bubble nucleation is initiated might be determined by further saturation measurements or by including nucleation into a transport model, cf.⁴³ Furthermore, the "Tafel-region" can be interpreted more accurately by taking the concentration increases due to current density dependent gas crossover into account.⁵⁷ Moreover, the breakdown results can be improved by approximating the mass transport overpotentials due to the (oversaturated) solute transport in catalyst layers, which strongly depend on their structural properties.⁵⁷ All these so far neglected effects may explain (part of) the discrepancy between breakdown and model results. Further work will shed light on it.

Representative volume elements were found for all investigated materials regarding their porosity¹³ and experimental investigations showed a Darcian, path-connected flow in similar PTLs.^32,33 Therefore, if adequately parametrized, the adopted continuum porous media model is expected to describe well the fluid states in PTLs at changing thermodynamic and flow conditions. This is reinforced by the agreement between model results and experimental data for the T10 material at different current densities and operating pressures, cf. Figs. 4 and 6. Non-continuum effects (e.g. limiting throats, as discussed above) can affect the transferability of parameters between materials and might have thus contributed to the fact that no match between model results and experimental data was obtained for T5 and T20 so far, yet previously discussed parameter uncertainties and inaccuracies in the breakdown method are expected to play an at least equally important role.

Let assume for the operating the conditions: i = 1 mA cm⁻², ϑ_op = 50 °C and p_op = p_ambient. Moreover, the following holds for the T10 material: transport path lengths $l\geqslant 1\,\mathrm{mm}$, porosity ϕ ≤ 0.37 and entry pressure ${p}_{e}\leqslant 5\,\mathrm{kPa}$.^6,13 The diffusion coefficient can then be estimated as

$$\begin{eqnarray*}\begin{array}{rcl}{D}_{\mathrm{pm},l} & \leqslant & {D}_{\mathrm{pm},l}(\phi =0.37,{S}_{l}=1)={\phi }^{4/3}\,{S}_{l}^{10/3}\,2.2\\ & & \times \ {10}^{-9}\,{{\rm{m}}}^{2}\,{{\rm{s}}}^{-1}\displaystyle \frac{{T}_{\mathrm{op}}}{298.15\ {\rm{K}}}\lt 6.4\times {10}^{-10}\,{{\rm{m}}}^{2}\,{{\rm{s}}}^{-1},\end{array}\end{eqnarray*}$$

see Eq. 15 and Ref. 44. The concentrations are calculated according to Henry's law 6. Their averaged (over Σ_a or Σ_c) gradients are

$$\begin{eqnarray*}\begin{array}{rcl}\parallel {\rm{\nabla }}{c}_{l}^{{{\rm{O}}}_{2}}\parallel & \leqslant & \displaystyle \frac{{c}_{l}^{{{\rm{O}}}_{2}}({p}_{g}={p}_{\mathrm{op}}+{p}_{e})-{c}_{l}^{{{\rm{O}}}_{2}}({p}_{g}={p}_{\mathrm{op}})}{l}\\ & \leqslant & \displaystyle \frac{0.853\,\mathrm{mol}\,{{\rm{m}}}^{-3}-0.806\,\mathrm{mol}\,{{\rm{m}}}^{-3}}{1\,\mathrm{mm}}=47\,\mathrm{mol}\,{{\rm{m}}}^{-4},\\ \parallel {\rm{\nabla }}{c}_{l}^{{{\rm{H}}}_{2}}\parallel & \leqslant & \displaystyle \frac{{c}_{l}^{{{\rm{H}}}_{2}}({p}_{g}={p}_{\mathrm{op}}+{p}_{e})-{c}_{l}^{{{\rm{H}}}_{2}}({p}_{g}={p}_{\mathrm{op}})}{l}\\ & \leqslant & \ \displaystyle \frac{0.660\,\mathrm{mol}\,{{\rm{m}}}^{-3}-0.623\,\mathrm{mol}\,{{\rm{m}}}^{-3}}{1\,\mathrm{mm}}=37\,\mathrm{mol}\,{{\rm{m}}}^{-4},\end{array}\end{eqnarray*}$$

with Eq. 7, ${\chi }_{l}^{\kappa }=\tfrac{{p}_{g}-{p}_{g,\mathrm{sat}}^{{{\rm{H}}}_{2}{\rm{O}}}}{{K}_{{\rm{H}}}^{\kappa }}$, ${c}_{l}^{{{\rm{H}}}_{2}{\rm{O}}}=54844.95\,\mathrm{mol}\,{{\rm{m}}}^{-3}$, ${p}_{g,\mathrm{sat}}^{{{\rm{H}}}_{2}{\rm{O}}}$ = $12351.27\,\mathrm{Pa}$, ${K}_{{\rm{H}}}^{{{\rm{O}}}_{2}}$ = $5.96087\times {10}^{9}\mathrm{Pa}$ and ${K}_{{\rm{H}}}^{{{\rm{H}}}_{2}}$ = $7.70828\times {10}^{9}\mathrm{Pa}$.^45,46 Finally, Eqs. 8, 10 and 13 give

$$\begin{eqnarray*}&&\begin{array}{l}\displaystyle \frac{\parallel {j}_{l,\mathrm{dif}}^{{{\rm{O}}}_{2}}\parallel }{\parallel {j}^{{{\rm{O}}}_{2}}(i=1\,\mathrm{mA}\,{\mathrm{cm}}^{-2})\parallel }\lt 0.006 \% \quad \ \mathrm{and}\\ \displaystyle \frac{\parallel {j}_{l,\mathrm{dif}}^{{{\rm{H}}}_{2}}\parallel }{\parallel {j}^{{{\rm{H}}}_{2}}(i=1\,\mathrm{mA}\,{\mathrm{cm}}^{-2})\parallel }\lt 0.046 \% .\end{array}\end{eqnarray*}$$