Understanding Impacts of Catalyst-Layer Thickness on Fuel-Cell Performance via Mathematical Modeling

In this article, a two-dimensional, multiphase, transient model is introduced and used to explore the impact of catalyst-layer thickness on performance. In particular, the tradeoffs between water production and removal through transport or evaporation are highlighted, with a focus on low-temperature performance. For the latter, a case study of an ultra-thin catalyst layer is undergone to explore how various material properties alter the steady-state and startup performance of a cell. The ﬁndings provide understanding and guidance to optimize fuel-cell performance with thin electrodes.

Polymer-electrolyte fuel cells (PEFCs) have emerged as a promising zero-emission technology for energy conversion due to their thermodynamic efficiency and high energy density. 1,2 However, to reduce cost, the amount of precious metal catalyst needs to be lowered. The most common strategy for such catalyst thrifting is to fabricate thinner catalyst layers. The prototypical example of this approach is the nanostructured thin-film (NSTF) electrode, which has several advantages compared to standard carbon-supported Pt electrodes. 2 These stateof-the art electrodes have demonstrated improved chemical stability, durability and desired specific power and activity, at the same time having high Pt mass activities. 2,3 Common to these and other lowloaded electrodes are the issues associated with water management in thinner electrodes. Typically, thinner electrodes are susceptible to severe flooding due to their inherently low water capacity and perhaps lack of hydrophobic zones. Such phenomena are particularly pronounced when PEFCs operate at lower temperatures or during startup.
Recently, several studies reported water-management mitigation strategies for thinner electrodes including modification of operating conditions and/or component morphologies to ensure successful startup and operation at low temperatures. For example, for NSTF electrodes, Steinbach et al. 4,5 reported a novel water-management scheme of increasing the pressure on the cathode to drive water to the anode; coupled with a different anode design, the limiting current density at low temperatures increased by a factor of four. Kongkanand et al. 6,7 demonstrated that water removal and storage capacity of the NSTF electrode can be significantly enhanced by placing an additional Pt/C layer between the electrode and microporous layer. The various empirical findings of NSTF as well as traditional supported thin electrodes 8 can be much better understood by examining the various tradeoffs engendered and complications arising from the use of thin catalyst layers.
Currently, the transport mechanisms of water removal behind possible mitigation strategies for thin electrodes are not well understood. In terms of modeling, both pore-scale and continuum models have been generated, although only the latter are germane to understanding cell water and thermal management. 9 Examples of continuum models include bilayer models, such as the one developed by Sinha et al. 10 with a membrane and water-covered portions inside of the electrode, although many of the parameters used are still not measured reproducibly. In terms of cell-level models with thin or variable-thickness catalyst layers, the most comprehensive were those of Nandy et al., 11 Balliet and Newman, 12 and Chan and Eikerling. 13 The former two focused on freeze and successful startup, which was also explored by others. [14][15][16] These works focused more on the key drivers of adding capacity in the membrane for water or getting water out of the catalyst layer before it freezes under isothermal conditions. These tradeoffs were not explored for higher temperature operation tradeoffs nor for more realistic startup conditions. In terms of Chan and Eikerling, they used a simplified model focused on water transport with ultra-thin-film catalyst layers and the related impacts of liquid-water transport, but did not examine what happens as a function of catalyst-layer thickness as the layer was treated as an interface.
In this article, we report on findings of water and thermal management as a function of catalyst-layer thickness, with a focus on the complications with thin electrodes, using a two-dimensional, transient, nonisothermal, full-cell model. It is the goal of this article to elucidate the key transport mechanisms and component design parameters to explain experimental trends and help guide future studies, as well as present a comprehensive modeling framework. The article is organized as follows. First, the model framework is presented with a discussion of all relevant physics and parameters. Next, the model is used to describe the impacts of catalyst-layer thickness on performance. Last, the ultrathin NSTF electrode is used as a case study to explore the impact of various material parameters in the hope of determining how one can increase its low-temperature performance.

Model Description
The PEFC model framework presented in this work is based on the model developed by Balliet and Newman 17,18 and previously by Weber and coworkers. [19][20][21][22] Numerous aspects of the model have    [2] Gas pressure Mass Fractions Membrane chemical potential Electric current density ∇ · i 1 = −i r xn h [10] Plate, GDL, MPL, CL Electric potential i 1 = −σ ef f ∇ 1 [11] Plate, GDL, MPL, CL Temperature α ε α ρ αĈ p,α [12] Plate, GDL, MPL, CL, M been modified and hence a full description of geometry and physics is presented here. The model is a two-phase, two-dimensional, cross-sectional sandwich model where model domains are shown in Figure 1. The modeling domains include membrane (PEM), cathode and anode plate (Plate), channels (CH) and lands (L), gas-diffusion layers (GDLs), microporous layers (MPLs), and catalyst layers (CLs). The reference point for the x-direction is set at the aCL|PEM interface; for the y-direction, the model spans from the mid-plane of the channel, where the reference is set to 0, to the mid-plane of the land. The base-case model subdomain dimensions and type of materials used are shown in Table I. The model accounts for different phases present in a particular domain as listed in Table I. For the study in this paper, the catalyst-layer thickness was varied; for the NSTF case study, it is set to the measured value of 0.5 μm. 2,10 Governing equations.-The governing equations, model boundary conditions and associated source terms are summarized in Tables  II through IV. Mass, electronic and ionic charge, and energy are conserved within the model with the use of second-order partial differential equations. The various parameters are given in Table V and  Table VI. Below, we discuss each major phenomenon in turn.
Energy balance.-The energy balance is solved on all of the domains with Eq. 12, where the model assumes local thermal equilibrium between phases to solve for temperature. On the right-side of Eq. 12, heat source terms account for water evaporation/condensation,Q v , joule heating, Q jle , and heat due to electrochemical reaction, Q r xn , which are tabulated in Table IV. In this work, experimental values are used extensively whenever possible as shown in Table V. The GDL (MRC 105 used in this modeling study) thermal conductivity in the in-plane direction is reported to vary with liquid-water saturation and is incorporated into the model as 50 k T,ef f,G DL = 0.187 + 0.524 (1 − exp (S L /0.697)) + 0.11 (1 − exp (−S L /0.0348)) [13] For the energy-balance boundary conditions, the temperatures are set at the outside edge of Plate as reported in Table III. The thermal contact resistance at the BP|GDL interface is incorporated within the model as [14] where, R G DL|B P = 2 × 10 −4 [m 2 K /W ] is the thermal contact resistance. 28 Gas convection and diffusion.-For the gas species within the porous domains (CLs, MPLs and GDLs), Darcy's Law (Eq. 3) and conservation of mass (Eq. 4) describe convection, whereas Stefan-Maxwell equations (Eq. 5) describe diffusion. The effective gas permeability, k ef f , of the porous domains is assumed to be a product of saturated permeability, k sat , and relative permeability, k r,G which depends on local liquid saturation value, S L k ef f = k sat k r,G = k sat (1 − S L ) β G [15] The exponents, β G , and saturated permeabilities for the porous domain are reported in Table V. The effective diffusion coefficient, D ef f i j , accounts for gas-species transport through porous and tortuous pathways of porous media and is described as 29

Variable Boundary condition Boundary
Liquid pressure, P L P L = P G for P L ≥ P G −n · N L = 0 for P L ≤ P G Anode CH|GDL Cathode CH|GDL where tortuosity, τ G is defined as ε −βτ G and porosity of gas-phase is related to that of porosity of domain, ε, through: ε G = (1− S L )ε Because Eq. 5 is an inverted form of Stefan-Maxwell equations, the binary diffusion coefficients,D i j , depend on mixture composition and are related to the diffusion coefficients used in Stefan-Maxwell equations, D i j , through the following expression for a three-component system 30 [17] where ω and y are the mass and molar fractions of gas species, respectively. The source terms in the conservation Equation 4 are due to evaporation/condensation in all of the porous domains and reactant depletion in the CLs as reported in Table IV. The gas pressure and mass fractions are specified at the anode and cathode CH|GDL interfaces. The mass fraction of oxygen is computed from the sum of all mass fractions equals one. In this paper, only fully saturated (RH = 100%) gas feeds are considered.
Liquid-water and multiphase aspects.-For liquid water within the porous domains (CLs, MPLs and GDLs), Darcy's Law (Eq. 1) and conservation of mass (Eq. 2) describe convection. The effective liquid permeability, k ef f , of the porous domains is treated in a similar manner to that of gas permeability [18] where S irr is irreducible saturation accounting for isolated regions of liquid water that do not participate in capillary transport. The exponents, β L , for every porous domain are reported in Table V. For numerical reasons, a negligible permeability at saturation equal to zero is added. Although the reported saturated permeability in GDLs is on an order of 10 −12 m 2 , 31 to realize a realistic pressure drop, a saturated liquid permeability of 10 −16 m 2 was used in the model. This order of magnitude was also adopted in a recent modeling work by Alink and Gerteisen. 32 The discrepancy arising from the oversimplification of the complex capillary-dominated transport engendered by using Darcy's law; 9,33,34 more complex models should be used as discussed recently. 35 Using this value the model's pressure drops through the GDLs are on the order of 1 to 10 kPa, consistent with the experimentally measured breakthrough pressures for different GDLs. 36,37 The liquid-water source terms are due to water production, water condensation/evaporation, and membrane sorption/desorption. The boundary condition at the CH|GDL for liquid pressure is not trivial and requires additional attention as we believe it to be of crucial importance. Previous approaches include setting the liquid pressure at this boundary to one that corresponds to a small saturation value (about 0.1) or to recast equations in terms of saturation; 38 however, such an approach is not physical as pressure is the thermodynamic driving force for liquid flow. In this manuscript, a liquid-pressure boundary condition is adopted using the approach of Weber and Newman 39 as Table III shows. For the base-case simulations, P L at the CH|GDL boundary is set equal to P G . The model assumes no reservoir of liquid water is present in the channel; therefore, the only physically meaningful direction of liquid-water flux is out from the GDL and into the Journal of The Electrochemical Society, 163 (7) F691-F703 (2016) If the direction of the flux of water is into the GDL, then a no-flux boundary condition is set. In several simulations reported herein, P L was set higher than P G at the CH|GDL boundary to simulate a GDL with a higher droplet adhesion force, 40 although the same allowance for only fluxes out of the GDL in the liquid phase was maintained. It is well known that droplets along this interface results in higher dynamic liquid pressures before the droplet is removed. [41][42][43][44][45][46] For multiphase properties, the saturation is determined using water-retention curves that plot the saturation as a function of the capillary pressure, In the model, a look-up table is created from the experimentally measured curves with liquid water imbibition/withdrawal technique 47,48 using a cubic-spline interpolation function. For the CL, measured data for a standard carbon-supported CL is used for all CLs including NSTF electrodes. 49,50 Figure S1 of the SI shows the water-retention curves used in the model. It should be noted that although water-retention curves are useful, there is some debate of their applicability when applied locally, 33 a topic of future study.
To relate the gas and liquid pressures, water source/sink terms are used as in Tables I and IV, where the membrane water sorption/desorption coefficient is water-activity dependent as reported by Kientiz et al. 27 and shown in Table VI. To determine the chemical po-tentials, thermodynamics is used. Following the derivation of Balliet, 17 and assuming ideal-gas behavior, constant composition of diluents in the gas phase, and constant heat capacities, the chemical potentials of vapor, μ V , and liquid water, μ L , can be expressed in terms of gas and liquid pressures and temperature relative to those at the triple point [21] whereH α,t andC α,L are molar entropy and specific heat capacity of water in phase α, respectively,V L is specific volume of liquid water, and y V is the molar fraction of vapor. One can set the two chemical potentials equal to each other to calculate the vapor pressure, which can be corrected for meniscus effects using the Kelvin equation 51 ln P sat where P sat,0 L is saturation vapor pressure when no diluents are present and P sat L is the actual saturation vapor pressure.  Vapor-equilibrated transport coefficient, Membrane.-In the membrane, transport of water and protons is considered. Eqs. 6 and 8 describe conservation of ionic charge and mass inside the membrane and CL domains, respectively. The driving forces for ion and water transport in the membrane are gradient of chemical potential of water, ∇μ 0 , and a gradient of ionic potential, ∇ 2 , as shown by Eqs. 7 and 9. A hybrid approach is used here, as proposed by Weber and Newman, 52 where the membrane has two transport modes: vapor-and liquid-equilibrated. A continuous transition is assumed between the two modes using a fraction of expanded channels, S. The membrane transport model has six transport parameters that are dependent on local conditions such as T, RH and water content, λ as listed in Table VI. Explicit membrane swelling effects are neglected in this model treatment since the feeds are fully humidified. The membrane used is a 24 μm thick membrane with equivalent weight of 850 g/mol developed by 3 M. The fit for a vapor-equilibrated water content, λ V , as a function of water activity, a 0 , was created using an in-house experimental data as reported in Table VI. The ionic conductivity of the membrane as a function of temperature was determined with fits from in-house measurements for 50% and 70% RH and the experimental data from Bosnjakovic et al. 53 for 90% was best fit by an activation energy of 14.5 kJ/mol, which is consistent with the activation energies reported for similar membranes. 54 It should be noted that the presented membrane model implicitly accounts for thermal-osmosis 55,56 through the temperature corrections to the chemical potentials (see Eq. 20 and 21). In this scheme, under the presence of temperature gradients, liquid water moves from the colder to hotter side of the membrane -in a direction that increases the entropy. 57 No-flux boundary conditions are applied on the outside surfaces of the CLs. In the model, it is assumed that water sorption/desorption happens inside the CL domain, where the water fluxes are defined as Electron transfer.-Conservation of charge applies to electron transport as shown by Eq. 10, where the electric current density, i 1 , is described by Eq. 11 using the electronic conductivities reported in Table V. The outside edge of the anode Plate is kept as a reference at 0 V, whereas either potential (potentiostatic mode of operation) or current (galvanostatic mode) was applied at the outside edge of the cathode Plate, depending on the operation mode.
Reaction kinetics.-For the reaction kinetics in the anode, a standard Butler-Volmer equation was used, [23] where, j H O R,0 , is the exchange current density, α H O R,a and α H O R,c are the anode and cathode symmetry coefficients, and the anode overpotential, η H O R , is defined as the difference between electronic, 1 , and ionic, 2 , potentials minus the standard half-reaction potential, The other relevant parameters are reported in Table VII.
For the cathode, the oxygen reduction reaction (ORR) is modeled with a double-trap kinetic model initially developed by Wang et al. 58,59 and built upon by the other recent studies. 60,61 The advantages of this kinetic model is that it predicts the observed doubling of the Tafel slope, 62,63,64 does not assume a reaction rate-determiningstep, and at the same time predicts oxide coverage values reported experimentally. 65 Depending on the operating conditions (potential), the ORR precedes through either or both adsorption pathways with four elementary reactions. The two pathways -O-adsorption and OHadsorption, as well as four elementary reactions (with labels) are [25] With the steady-state approximation, the kinetic current for the ORR is expressed as a combination of elementary reactions (RT + RA, RA + DA), or as a single RD-step where, j * is a reference prefactor, G * R D and G * −R D are the potential-dependent activation free energies of forward and backward reactions, respectively. These activation energies can be described in terms of activation, G 0 R D , and adsorption, G 0 O H , free energies: Agglomerate and thin-film model.-To capture the heterogeneous structure of the CL, an agglomerate model is used, which is only implemented on the cathode side, as mass-transport losses on the anode side due to hydrogen dissolution into agglomerates are assumed to be negligible. Agglomerates are assumed to be of constant size and covered by a thin continuous ionomer film. The agglomerate current is given by 66 where r agg and δ agg are the agglomerate radius and thin-film thickness, respectively. P O2 /H O2,N is the ratio of oxygen pressure and Henry's constant at the thin-film|pore interface, which is also the oxygen concentration,V agg = r 3 agg /(r agg + δ agg ) 3 is the ratio of volume of agglomerate core to the volume of entire agglomerate. The agglomerate current is related to the volumetric catalyst layer current, i r xn . Assuming a first-order oxygen dependence, the analytical solution for the reaction effectiveness factor, E r [28] where, φ L is Thiele's modulus for chemical reaction: [29] The reaction rate at the surface of agglomerate is: where (1 − ε C L )V agg is the active area scaling factor. Numerical method.-The above governing relationships were solved using COMSOL Multiphysics 5.0/5.1 (COMSOL, Inc., Palo Alto, CA) coupled to Matlab R2014a with steady-state and transient solvers. The mesh was generated with 5,388 triangular and 2200 quadrilateral elements, where increased mesh density was introduced within and near CL domains. The fully-coupled physics were solved with MUMPS general solver, with 110,563 degrees of freedom. Parametric sweeps were run through a coupled Comsol/Matlab interface, where defined variables and functions were passed between the two programs. This enabled ease of data storage in a form of structures through a Matlab interface.

Results and Discussion
Model validation and polarization behavior.-Since the focus of the study is on the impact of catalyst-layer thickness, especially thin ones, the model was first validated by comparison to data from NSTF cells. In the following sections the model results and relevant discussion are presented. Figure 2 shows model prediction for experimental data for NSTF PEFCs at two different temperatures. The experimental data is for two different cell builds with Baseline GDLs and operating conditions listed in Table VIII. The model shows a good agreement with the polarization curves, predicting high current densities, 2 A/cm 2 at 0.45 V at 80 • C and almost an order of magnitude reduction of limiting current, 0.21 A/cm 2 at 0.45 V for PEFC operating at 40 • C. This temperature sensitivity is attributed to poor water management and is discussed in more detail in the later sections and is a key factor that the model should capture. It should be noted that the high-temperature polarization curve for NSTF is also representative of thicker, traditional Pt/C catalyst layers.
Furthermore, it is important for the model to be able to predict startup of the PEFC at low temperatures, as this could be greatly impacted by catalyst-layer thickness or water capacity. The transient simulations were run in galvanostatic mode, where a constant current of 0.25 A/cm 2 was applied and cell voltage monitored over 400 s. Figure 3a shows start-up transients for two types of cell builds -adiabatic  and isothermal. Adiabatic cell architecture developed by United Technologies Corporation (Hartford, CT) simulates the cell inside a PEFC stack, 69 where major heat sinks due to the endplates do not exist as they do for a typical single-cell or isothermal test. For the adiabatic case, the model predicts a linear increase in cell temperature over time from 40 • C initially to 51 • C at time of 400 s (Figure 3b), whereas the experimental data show increase only to 46 • C; the discrepancy arises from the fact that the model uses a truly adiabatic or no-flux boundary condition for energy whereas experimentally some losses remain. Moreover, in the experimental set-up the thermocouples were placed in the locations under the gas-flow in the channel, thus possibly under-predicting the cell's temperature due to gas convective cooling. Comparison of the model with experimental data shows good agreement, where, during the initial adiabatic cell heating, the potential drops but recovers within the initial 100 s, after which the potential steadily increases and startup of the cell is successful. Whereas for the isothermal cell, the potential continues decreasing, although at a less steep rate, until cell shutdown happens at 450 s. This initial voltage decay is due to a time-lag in membrane humidification and low cell temperatures. The potential recovery for the adiabatic cell is due to higher temperature phenomena including property changes occurring. Figure 2, there is quite a performance discrepancy with temperature for the very thin NSTF electrodes. However, one would expect that performance difference to change as a function of catalyst-layer thickness. To study this issue, simulations are run at various operating temperatures and catalyst-layer thicknesses as shown in Figure 4. In these simulations, the specific reactive surface area is held constant and thus the amount of catalyst varies with catalyst-layer thickness, similar to what occurs in reality. As shown in Figure 4, decreasing the loading plays a much less significant role in terms of affecting performance compared to the difference in performance due to temperature change. The results demonstrate several interesting features. First, at high temperature, the impact of thickness is minimal with the major effect being the decreased amount of catalyst loading, although this is not substantial. For the very thin catalyst layers, flooding occurs which limits the performance. It is interesting that there appears to be a slight peak around 5 micrometers due to minimizing some of the mass-transport issues within the catalyst layer and especially the impact of having a higher average temperature in the catalyst layer due to the higher volumetric reaction rates as discussed below. This result also shows that under these conditions only about 5 micrometers is well utilized. Figure 4 also demonstrates that for thicker catalyst layers, the dependence on temperature is not that strong, due to the ability of catalyst layers to minimize the impact of saturation. In fact, only for thicknesses be-low 5 micrometers is there a significant change in the temperature response, and for all thicknesses, that change occurs around 50 • C. The reason for this change is that this temperature represents when the cell moves from mainly liquid to vapor water transport through phase-change-induced (PCI) flow. 22,71,72 Thus, at this temperature, the combination of thermal gradients and vapor pressure result in water evaporating in the catalyst layer and moving along the vapor-pressure gradient and condenses in the channel or at the land. 71 Such impact is shown in the higher water saturation at the land in Figure 5a. It should be noted that due to this change, the intermediate temperature range around 50 • C also results in the hardest numerical convergence for the model.

Impact of catalyst-layer thickness.-As shown in
The above issues are demonstrated in Figure 5b, where it is clear that flooding is not significant at higher temperatures due to PCI flow. At lower temperatures, flooding occurs but the thicker catalyst layers can accommodate this without impacting performance significantly. Figure 5a shows that the liquid saturation at the anode is quite high, but due to the pure hydrogen feed and fast diffusion, this does not impact the performance. The very high liquid saturation is driven by the fact that there is a substantial amount of water being transferred from cathode due to the thinness of the catalyst layer localizing the water generation and thus liquid pressure within the smaller volume. This is seen in that the net water flux, β, is slightly negative (−0.003) for 0.5 μm and 40 • C, and increases with thickness to −0.0005 for 10 μm as shown in Figure 5c. The direction of the net membrane water flux, N mem w , is defined as positive when water is transported to the cathode and negative when it is transported to the anode.
Beyond the impact of thickness on steady-state performance is its influence on startup, where there is a competition between water production and flooding and water removal and heat production. To explore these interactions, transient simulations were conducted. Previous work has highlighted that these tradeoffs are key for understanding operation and successful startup, especially when time-dependent phenomena occur such as phase-change effects including evaporation and freezing. 11,12,[14][15][16]73 Figure 6 shows simulation results of adiabatic startups for various catalyst-layer thicknesses; it should be noted that the timescales are sensitive to input parameters including heat capacities, but the trends remain comparatively the same. As seen in the figure, thicker catalyst layers result in better performance due to their lower saturation caused by the lower volumetric generation rate. This lower rate also results in a faster approach to stable performance, although not necessarily a faster heat up, and subsequent performance increase as seen (see Figure S3 in the SI) in the slightly smaller slope increase with the thicker catalyst layers. For thin catalyst layers, a higher volumetric rate causes a deeper and significantly longer performance decrease while the tradeoffs between water and heat production and water removal balance. Thinner layers also exhibit a continuous decrease in performance when simulated isothermally at 40 • C (see Figure 10) and a lower sensitivity to thickness in agreement with Figure 4.

Property sensitivities for thin catalyst layers: NSTF case study.-
As mentioned, NSTF is the prototypical thin-film catalyst layer with about 0.5 μm thickness, and thus serves as a good test case for exploring the sensitivities of various properties on low-temperature performance using thin catalyst layers. As shown in Figure 2 and elsewhere (see also Figure 4), these electrodes demonstrate extreme sensitivity to temperature, which can be mitigated by moving water from the cathode to the anode and thus relieve cathode flooding. 4,5 In this section, we explore several properties that could impact such a mechanism and explore their relevance. It should be noted that the focus here is on material properties rather than operating conditions, since these are assumed to be specified by the system and stack and thus not easily altered. For example, the most obvious choice to force water through the anode is to increase the cathode versus anode pressure (or even go to vacuum pressures), 5,74 but which might not be practical at the system level.
Since the desire is to move the water out of the anode, the properties explored are GDL permeability and interfacial liquid pressure, the membrane sorption rate and diffusivity, and the electrode hydrophobicity. While cathode properties and changes (e.g., thicker or more resistive MPLs) could likewise promote more water out of the anode, 45,75-77 these often result in additional cathode-catalyst-layer flooding (which is not the case for the increased cathode gas pressure mentioned above) and so are not considered. A tornado plot showing the sensitivities of the various parameters is shown in Figure 7. As seen in the figure, the most sensitive parameters are the interfacial   Tables V and VI. anode GDL|channel liquid pressure and CL liquid pressure, whereas some of the others have to be varied orders of magnitude in order to impact performance. The baseline current for the baseline properties used in the model is 0.95 A/cm 2 (at 0.4 V), and the current decreased or increased depending on each individual change in property value. For example, reducing a property by 10 (×0.1 in Figure 7) results in currents of 0.91, 0.5, and 0.75 A/cm 2 for anode liquid permeability, membrane sorption rate, and membrane diffusivity, respectively. Of course, different material properties are obtainable with different degrees of ease. From this comparison, the membrane sorption rate is the most sensitive parameter. As discussed below, the increased performance tracks strongly with an increase in the net water flux being driven out the anode than the cathode (i.e., more negative β). Each parameter is discussed below in turn.
Anode-GDL liquid pressure.-The liquid pressure at the CH|GDL interface arises from the existence of discrete water droplets that must be removed by the flowfield either by contact with the lands or by detachment due to gas flow in the channel. This droplet detachment is strongly coupled to the adhesion force of the material on the GDL surface. 40,41,46,78,79 Incorporation of this physics into continuum cell models remains a challenge. 9 However, such effects should essentially impact the water liquid pressure at the boundary, which is the pressure of the droplet at detachment, assuming Darcy's law remains valid. Furthermore, it is suspected that this value is dependent on GDL type and perhaps even history. 78 To explore this impact, the liquidpressure boundary condition was varied, with the impact on current shown in Figure 7 (see also Figure S4 of the SI). As the pressure is decreased from 160 to 100 kPa, the cell current density increases 5× to almost its higher-temperature performance (see Figure 4) due to a much lower amount of flooding within the cathode (as shown in Figure 8a (compare to Figure 5a)) caused by the lower liquid pressures resulting from the increased flux to the anode and lower anodeboundary liquid pressure. This net flux of water is shown in Figure  8b, where β (see Eq. 31) is plotted as a function of boundary liquid pressure.
It is apparent that the best performance occurs with the most negative β, corresponding to the most water out of the anode. For a comparison of the actual water fluxes, see Figure S5 in the SI, although it should be noted that these simulations are at 0.4 V, and thus the net water flux varies for different liquid pressures, which is one reason why β is a good comparison parameter since it normalizes this effect. At higher liquid pressures, β steeply increases, reversing sign at 161 kPa, indicating net water transport from anode to cathode. At 180 kPa, β is at its limiting value, corresponding to a zero liquid flux as shown in Figure S5 in the SI. For smaller anode-boundary liquid pressures, β asymptotes to a value of −0.01, and further decreasing anode liquid pressure will not enhance water removal through the anode side. Under these conditions, droplet removal into the channel is no longer a limiting step in anode water removal but some other  mechanism (e.g., water transport through anode GDL or membrane) becomes more important.
From a design consideration, a GDL with a low adhesion force is preferred, where liquid-droplets can be detached with lower resistance. Several design mitigation parameters to enable this include GDL or land surface treatments, GDL architectures with morphological changes including perhaps modulated porosity, 80 cell compression level, higher gas flowrates, or different channel/land dimensions. These modifications will not necessarily benefit GDL-MPL or GDL-CL interfaces, where low water breakthrough pressures are desirable for water removal from the CL. To accommodate GDL-CH and GDL-MPL/CL interfaces, a spatially inhomogeneous structure is perhaps preferred, where the GDL properties are modified chemically and structurally to aid in water removal from the CL while simultaneously having a small droplet adhesion force on the other side of the GDL.
Anode-GDL permeability.-The permeability of the GDL is tightly coupled to the GDL structure and surface treatment. GDL hydrophobization can possibly reduce the porosity and permeability of a GDL, as reported in the literature. [81][82][83][84] In addition, compression and landchannel effects can have similar effects on transport parameters. 37,[85][86][87][88][89][90][91][92] Locations of GDLs under land are more compressed, thus having lower porosity and hence lower permeability. To explore the impact of permeability, the anode liquid saturated permeability was decreased by factors between 10 and 100, and increased by factors 10 and 100 with the results as shown in Figure 7 and Figure 9a. Increasing permeability by a factor of 100 did not show dramatic effect on current density, as it increased to 1.03 A/cm 2 , indicating that the anode GDL permeability is not a performance-limiting factor. As the permeability decreased by two orders of magnitude from the baseline, the current decreases from 0.95 to 0.63 A/cm 2 . This 37% reduction in current density is significant; however, even with high compressions or extreme PTFE loadings (>40 wt-%) a 100× decrease is not expected. A decrease in one order of magnitude, which is physically justifiable, results in only a minor 9.4% reduction in current density. The change in current density tracks that of β as shown by Figure 9a as it did for the liquid-pressure case above (see Figure 8). The 100× decrease in permeability results in increased liquid pressure such that the anode pressure drop increases across the GDL (up to 30 kPa), significantly changing the saturation throughout the cell (Figure S6c of SI) and liquid water flux ( Figure S6d of SI), with a significant impact on anode flooding.
Membrane sorption and transport coefficient.-A key aspect of removing water through the anode is getting the water across the membrane. In terms of transport properties, the sensitivity is not that strong as shown in Figure 7, which is consistent with both the thinness of the membrane (24 μm) as well as the fact that transport of water through the membrane is seemingly interfacially dominated. For example, reducing the transport coefficient by an order of magnitude decreases the current density to 0.74 A/cm 2 , whereas decreasing the membrane sorption rate by an order of magnitude results in a current density decrease to 0.5 A/cm 2 . In terms of percentage change, these constitute 22 and 47% decreases, respectively. Increasing the transport coefficient does not affect the current density, since it is no longer limiting (see Figure S8 in SI). This is especially true for modern cells that utilize membranes on the order of 10 μm, where the bulk property is no longer that significant; however, interfacial effects could still limit water transport. It has been shown that the interfacial mass-transfer coefficient for water transport can be significant relative to the transport coefficient, especially for drier membranes. 27,93-100 Only a few models assume non-equilibrium and incorporate interfacial resistance for species transport. [101][102][103] To explore this effect in more detail, the baseline value adopted from Kientiz et al. 27 and reported in Table  VI was reduced by factors of 10 (×0.1) and 20 (×0.05) as shown in Figure 7 (and in Figure S8 of the SI). As the figures show, as the water sorption coefficient decreases, the current density decreases too. Further reducing the sorption rate has no significant impact on the current density because for these low values, most of the water produced is removed through the cathode side as evidenced by the dramatic change in β and the water fluxes shown by Figure 9b. Increasing sorption coefficient by factor of 10 increases current density from baseline 0.95 to 1.31 A/cm 2 , further insignificant increase in current is observed for a factor of 100 as shown by Figure 7.
In addition, one expects that the membrane sorption coefficient will impact the transient response as well. This is seen in Figure 10, where the two order-of-magnitude decrease in sorption rate results in a dramatic cell shutdown due to flooding since all of the produced water remains in the catalyst layer immediately upon startup. Particularly, Figure 10 shows simulations for 0.5 or 5 μm catalyst layers at a current density of 0.25 A/cm 2 , where the temperature outside of the bipolar plate is set to a constant 40 • C. For the 0.5 μm thick electrode, the potential drops to 0.45 V after the first 10 s and then continues decreasing to 0.42 V at 200 s, when the simulation was terminated. For the 5 μm thick catalyst layer, the potential recovery occurs faster (5 s) and a higher cell potential is maintained (0.49 V). Moreover, the potential decay is less steep, amounting only to 0.01 V for 200 s of simulation time. Such limited water capacity of thin catalyst layers means that startup, especially from low and subzero temperatures, requires membranes that have water capacity and can uptake water quickly as mentioned above.
Catalyst-layer hydrophobicity.-In almost all of the above simulations, cell performance became limited at the lower temperatures due to cathode catalyst-layer flooding. Thus, if one can alter the hydrophobicity of the catalyst layer, one could minimize the onset of flooding and thus increase performance. To examine this, simulations were run with different water-retention curves by shifting them to the left (positive shift) or right (negative shift) by a specified capillarypressure value (see Figure 11) for steady-state (see Figure S9 of SI) and transient conditions. As shown in Figure 7, this results in a significant change in the current density due to the varied levels of liquidwater saturation inside the CL. Increasing hydrophobicity of the CL (shifting water-retention curve by −20 kPa) results in current density increase to 1.11 A/cm 2 , which is 17% higher than the baseline. On the other hand, for highly hydrophilic CL (shift by +20 kPa) the observed current was 0.66 A/cm 2 , which is 31% decrease in current. For these highly hydrophilic CLs both anode and cathode are flooded (see Figure S9c).
Furthermore, as shown in Figure 11a, there is also a significant impact on the startup response of the cell. As the hydrophobicity increases, the cell's potential recovers to a higher value after an initial 10 s decay. The improvement in the potential recovery for the more hydrophobic layers is mainly due to their lower saturation, as shown by Figure 11b. Compared to a baseline case, the saturation decreases dramatically for a 30 kPa shift (from 0.92 to 0.58 at simulation times larger than 10s). However, the increased hydrophobicity also results in a more positive β (see figure S10 of the SI), thus more water out the cathode than the anode due to changing the driving force of water through the membrane. It should be noted that there is still debate as to whether retention curves or water-film formation are more limiting in CLs. Furthermore, it should be noted that it is not trivial to change the water-retention curve by 30 kPa, seeing as GDLs are typically in the range of 10's of kPa before complete saturation, and even MPLs are in the range of 100 kPa or so, although this is still an area of research and depends on fabrication method, cracks, etc. 31,47,49,50,[104][105][106][107] Compounded on this is in the need for catalysts at sufficient loadings in the thinner electrodes, thereby resulting in extended hydrophilic surfaces. However, hydrophobicity changes could explain recent observations of the impact of interlayers with NSTF showing increased performance.

Summary
In this article, the startup and steady-state performance of polymerelectrolyte fuel cells with varying catalyst-layer thickness was examined using a two-dimensional, two-phase mathematical model that incorporates the relevant physics. The model was validated with both steady-state and transient startup data using thin-film catalyst layers based on the nanostructured thin-film architecture. It was found that for catalyst layers thinner than about 5 μm, there is a significant dropoff in performance with temperatures below 50 • C, where phase-change-induced flow can no longer remove a majority of the produced water in the vapor phase. To mitigate the flooding of these thinner electrodes, various material properties were explored. The key finding was the need to remove more of the water through the anode than the cathode, thereby decreasing the cathode saturation. Such a mechanism can be enhanced by decreasing the anode gasdiffusion-layer / gas-channel boundary pressure for the liquid phase, assumed to be caused by water-droplet interactions there. Similarly, changes to the cathode catalyst-layer water-retention behavior can avoid flooding, although enacting such changes experimentally is not trivial. Less sensitive were changes to gas-diffusion-layer and membrane permeability. Of middle sensitivity was the membrane sorption coefficient. Overall, the findings and analysis presented help to explain the impacts of changing catalyst-layer thickness, especially at Figure 11. Transient nonisothermal simulation of a) isothermal startup for four different hydrophilicity/hydrophobicity levels in the cathode catalyst layer. The inset shows the associated water-retention curves.