Through-the-Membrane Transient Phenomena in PEM Fuel Cells: A Modeling Study

This paper presents a 2D, fully coupled and comprehensive transient model that accounts for micro-structural features of various cell layers. The model beneﬁts from state of the art sub-models for reaction kinetics and incorporates the polymer relaxation dynamics. Furthermore, a mixed wettability model is utilized to simulate the transient two phase conditions in the porous layers. The model is validated with transient experimental data under various conditions. A comprehensive simulation study is presented to investigate the impact of operating temperature and relative humidity on the transient response. The effects of cathode Pt loading and operation mode, i.e., current control versus voltage control, are also studied. The cell response is found to be dominated by water transport through its thickness. Additionally, it is found that reducing the Pt loading can inﬂuence the performance by changing the water balance in the cell, which has rarely been highlighted in the literature. In particular, at low temperature more water is transported toward the anode when the cathode Pt loading is reduced, since the resistance to water back diffusion is lowered with reduced thickness of the cathode catalyst layer. This trend is reversed at a higher temperature due to increased volumetric heat generation with reduced thickness. The model can help in understanding various transport phenomena and is expected to be useful for inspecting spatio-temporal temperature,

Some of the significant technological challenges to commercialization of polymer electrolyte membrane (PEM) fuel cells have been addressed in the past two decades through extensive theoretical and experimental research. As a result, the PEM fuel cell technology has advanced considerably. In particular, the amount of precious metal used in the catalyst layers (CL) has been significantly reduced while achieving remarkable performance improvements. 1,2 Nevertheless, important challenges related to cost and durability remain. While further material development is indeed beneficial, some of the existing issues may be addressed through effective control and hybridization of these systems. This will require a good understanding of the processes that govern the fuel cell dynamics. Moreover, in-depth studies of the transients can improve the current understanding of various electrochemical and transport phenomena. Therefore, there is a need for better understanding of the transient response to further enhance performance and lifetime of PEM fuel cells.
Arguably, the dynamic studies of PEM fuel cells have been overshadowed by the significant efforts dedicated to steady state measurements and modeling. However, fuel cell transient response has attracted some attention lately, as it can be used to elucidate and deconvolve complex transport phenomena. [3][4][5][6][7][8] Several transient models have also been proposed in the literature. [9][10][11][12][13][14][15] These models usually use simplified reaction kinetics and do not account for the microstructure of various cell layers and the anisotropic material properties. Therefore, these models typically do not have the required fidelity to allow detailed investigation of the transient phenomena affecting the cell performance. Accordingly, it is the main objective of this work to develop a transient model that captures the most salient features of the cell's dynamics across its thickness. Furthermore, we execute the model under a variety of operating conditions to delineate the critical transient phenomena that determine the overall cell dynamics. It should be noted that the model presented here is only helpful in understanding transient phenomena through the thickness of a unit cell * Electrochemical Society Student Member. z E-mail: jchen186@ford.com; tersal@umich.edu with small active area. In other words, the compressor and channel flow dynamics, along the channel redistributions, and stack thermal dynamics are not discussed here. Nevertheless, the results of this work may be extended along the flow channel and across multiple cells to study the transients at those scales. Specific to our objective is creating a model that, as much as possible, directly translates the physically measurable parameters and operating conditions into the knowledge about spatio-temporal distributions of critical variables such as temperature and water concentration in different layers. To this end, it is imperative to effectively capture the physical characteristics of the porous layers including the catalyst layers and diffusion media (DM). It is worth pointing out that, as has recently been shown, a representative elementary volume (REV) cannot be clearly defined in the through-plane direction for commercial DM 16 and the REV for the in-plane direction is on par with the land-channel sizes used in fuel cells. 16,17 Therefore, while the macroscopic models can capture the aggregate behavior, their predictions may significantly deviate from the local predictions by microscopic models. 16,18 This result bears significance, as it points to the need for more elaborate description of transport in the porous layers. Nevertheless, the excessively high cost of such simulations limits their application to very few conditions and a limited material set. Therefore, the macroscopic models remain the main tool to investigate the internal distributions in a full cell model and including some level of description of the microstructure is the approach adopted in this work to improve their prediction capability. This is achieved here by using a recently developed mixed wettability model for the porous layers. 19,20 To the best of the authors' knowledge, this is the first time that a full implementation of the mixed wettability model is being used in a 2D transient model. In addition to the mixed wettability model for porous layers, the presented model accounts for ionomer relaxation dynamics and CL micro-structure, which have been neglected in most of the previous models. Moreover, the effective material and transport properties of the different layers are identified through an extensive literature review on commercially available materials. Therefore, this model is expected to offer higher fidelity than the state of the art models for the purpose of studying transient phenomena that impact performance and durability. 21 The rest of the paper is organized as follows: first, the model formulation is presented along with a detailed review of the literature relevant to each sub-model to justify the choices made during model development. Simulation results and discussions are provided next, followed by a brief summary and concluding remarks. Model validation with experimental data from the literature is presented in the accompanying Supplementary Information.

Model Formulation
The modeling domain of interest is shown in Fig. 1. The model draws from prior work by Balliet et al., 10 Zenyuk et al., 22 and Zhou et al. 23 with modifications to several sub-models. Therefore, the model is presented here in its entirety for completeness. A complete list of model parameters is provided in the accompanying Supplementary Information.
Governing equations.-The comprehensive 2D model solves the following governing equations, where the various source terms are given in Tables I and II: ∂ (ρ l ε l ) ∂t = ∇ · ρ l K eff l μ l ∇ p l + S liquid , [2] ∂ (ρ g ε g ) ∂t = ∇ · ρ g K eff g μ g ∇ p g + S gas , [3] ε ion ρ ion EW ∂λ ∂t = ∇ · (N w,mb ) + S λ , [4] α ε α ρ α c p,α ∂T ∂t + u α · ∇T = ∇ · (k eff T ∇T ) + S T , [5] ∇ · (σ eff 1 ∇φ 1 ) = S e − , [6] ∇ · (σ eff 2 ∇φ 2 ) = S H + . [7] The first equation describes species transport in the microporous layers (MPLs), gas diffusion layers (GDLs), and CLs, where ε g is the porosity of the layer available for gas transport (ε g = ε(1 − s), with s being the liquid saturation and ε the compressed layer porosity), c g is the total gas concentration, u g = − K eff g μg ∇ p g is the velocity of the gas phase, and c i , D eff i , and x i denote the concentration, the effective diffusivity, and the molar fraction of species i, respectively. The first and second terms on the right hand side model the diffusive and convective fluxes, respectively, while the last term (S i ) is the relevant source term for the specific gas species (see Table I). In the anode, the equation is solved for water vapor molar fraction (x H 2 O ) and the hydrogen molar fraction (x H 2 ) is calculated by x H 2 = 1 − x H 2 O . On the cathode side, this equation is solved for water vapor (x H 2 O ) and oxygen (x O 2 ) and the nitrogen molar fraction is found from Equations 2 and 3 describe the mass conservation of the liquid and gas phases, respectively. These equations model the pressure drop of each phase in the CLs, MPLs, and GDLs according to Darcy's law. Here, ρ α , ε α , K eff α , μ α , and p α denote the density, volume fraction, effective permeability (relative permeability times absolute permeability), viscosity, and pressure of phase α. Note that ε g = ε(1 − s) and ε l = εs. Finally, S liquid/gas denotes the source term for the corresponding phase (see Table I). Note that liquid saturation that appears in these equations through the volume fractions is a variable that depends on capillary pressure. Therefore, closure equations that relate the saturation level to the capillary pressure are required. In this work, the mixed wettability model is used to derive water retention curves (liquid saturation vs. capillary pressure) as well as effective property values such as gas and liquid permeability for the different layers (see Mixed wettability model for porous layers section). Equation 4 governs water transport in the ionomer phase throughout the catalyst coated membrane (CCM). Therefore, its domain of application is the anode and cathode catalyst layers and the membrane. In this equation, ε ion , ρ ion , and EW denote the ionomer volume fraction, density, and equivalent weight, respectively, while λ is the dimensionless number that quantifies the water content in the ionomer, i.e., the number of water molecules per sulfonic acid group. Finally, S λ is the source term (see Table I) and N w,mb is the water flux in the ionomer phase across the CCM, which includes the effects of electro-osmotic drag (EOD), diffusion, and thermo-osmosis and is Table I. Mass conservation source terms.
Journal of The Electrochemical Society, 166 (7) F3154-F3179 (2019) calculated as follows: [8] where n d is the EOD coefficient, φ 2 is the ionic potential, F is the Faraday's constant, σ eff 2 is the effective conductivity in the ionic phase, D eff w,mb is the effective water diffusion coefficient in the membrane and D T ,mb is the thermal water diffusivity in the membrane. Note that thermo-osmosis is shown to drive water from the cold to the hot side for a hydrophilic membrane. 24 As a convention, a positive flux denotes water flux toward the cathode. The membrane water transport properties are given in Table III.  Equation 5 is the energy conservation equation, which governs the temperature distribution. In this equation, ρ α , ε α , c p,α , and u α = − K eff α μα ∇ p α are the density, volume fraction, specific heat capacity, and velocity of phase α, where α can be the gas, liquid, or solid phase. In addition, k eff T is the effective thermal conductivity (see Calculation of effective properties section) and S T denotes the heat source term (see Table II). Note that this equation captures heat transfer by conduction (first term on the right hand side) as well as convection (second term on the left hand side).
Equations 6 and 7 are the Ohm's law for electronic (φ 1 ) and ionic (φ 2 ) potentials, respectively. Here, σ eff 1 and σ eff 2 denote the effective conductivities of the respective phases and S H + /e − is the relevant source term (see Table II).
The source terms for mass conservation equations (Equations 1-4) are given in Table I. Here, M i is the molar mass of species i, i HOR/ORR is the volumetric HOR/ORR reaction current density, and S pc is the source term due to phase change and is given by where k evp/cnd denotes the rate of evaporation/condensation, a LG is the interfacial area between the liquid and gas phases (calculated by the mixed wettability model), p v is the vapor pressure, and p K sat (p c , T ) is the corrected saturation pressure that takes the Kelvin effect into account. The corrected saturation pressure is given by: [10] In the above equation, p c is the capillary pressure (p c = p l − p g ), R is the universal gas constant, ρ v is the density of water vapor, and p sat (T ) is the saturation pressure as a function of temperature given by: , [11] where T is in Celsius and the calculated pressure is in kPa.
In Table I, S ad denotes the source term due to water exchange (absorption/desorption) between the ionomer phase and the pore space and is given by: [12] where k ad is the interfacial water transfer coefficient (see Table III), δ CL denotes the CL thickness, and λ * is the dynamic variable for equilibrium water content in the ionomer (see Ionomer water uptake section). Note that water production with electrochemical reaction contributes to S λ . In other words, the produced water is assumed to be in absorbed phase. This is in agreement with the assumed structure for the CL in this work and has also been used by others. 9 The source terms for energy and charge conservation (Equations 5-7) are given in Table II. In the table, i 1 and i 2 are the electronic and ionic current densities, respectively: Moreover, H HOR/ORR denotes the reversible and irreversible heat of reaction given by: where HOR/ORR is the Peltier coefficient for HOR/ORR, 30 [17] where T is in Kelvin. Finally, H ad is the heat of sorption (due to water exchange between the ionomer and the pore space, i.e., water vapor) and is given by: 31 Ionomer water uptake.-It is imperative for any transient model of a PEM fuel cell to properly capture the dynamics of water sorption, desorption, and transport across the membrane. Historically,  29 and hydraulic 32 type models have been used for this purpose. However, there is abundant evidence in the literature suggesting that interfacial transport phenomena as well as swelling of the polymer backbone may play a significant role in water uptake and transport dynamics. 33,34 In particular, the gravimetric water uptake experiments conducted by Satterfield et al. have shown very long time constants for membrane water sorption while the desorption time constants were found to be an order of magnitude smaller. 34 They suggested that the sorption behavior may be explained through the contributions of interfacial water transport and stress relaxation in the polymer, whereas the desorption dynamics are dominated by the interfacial phenomena. Their experiments included a step in the humidity from fully dry to fully saturated conditions that resulted in significant relaxation behavior. Other studies have found much less pronounced impact of the relaxation dynamics when the membrane was subjected to smaller changes in the humidity conditions. 35 Similar results have been reported for ionomer thin films. 36 Dynamic vapor sorption (DVS) experiments by Kusoglu et al. have also shown relatively long water uptake times with a time constant that increases with membrane hydration. 37,38 Their results also indicate that the asymmetry between sorption and desorption is not as pronounced as that observed in Satterfield et al.'s experiments. Such significant difference was also challenged by Silverman et al. who found the desorption to be only about twice faster than sorption. 39 In situ measurements of membrane swelling and hydration by GM researchers have also found the hydration and dehydration rates to be similar. 40 Based on the preceding discussion, it stands to reason to incorporate the slow ionomer water uptake process into the model. Silverman et al. have developed a coupled transport and mechanical model that captures such phenomena. 31,39 However, adding the mechanical model will result in additional complexity that must be avoided for the purpose of performance modeling. Therefore, we simply use a dynamic variable to represent the quasi-equilibrium water content: 15,41,42 λ * = (1 − ϕ)λ eq + s relax , [19] where ϕ determines the contribution of relaxation phenomena to the ionomer water uptake (a value of 0.15 is used for the simulations in this work), and s relax is a variable accounting for the dynamics of stress relaxation. In particular, its dynamics are assumed to be first order: [20] where λ eq is the equilibrium water content given in Table III and τ is the relaxation time constant defined as: τ = exp(2 + 0.2λ). [21] Note that the time constant was chosen to vary with the ionomer water content in accordance with evidence in the literature for this dependence. 37 Moreover, the stress relaxation is supposedly a thermally activated process. 34,43 Therefore, it would be reasonable to assume an Arrehnius type temperature dependence for the associated time constant, which is not included here. Moreover, the effects of compressive stresses on membrane water uptake, 40,44 and the contentiously debated discrepancy between water uptake by ionomer thin film and bulk ionomer 36,37 are not taken into account in the model. Future parametric studies should aim at investigating these effects along with the effects of changes to the relaxation model parameters to understand their impact on the overall dynamic response of the cell.
Catalyst layer model.-Conventional catalyst layers of PEM fuel cells consist of Pt catalyst particles dispersed on carbon primary particles and (partially) covered by ionomer. The Pt particle size is in the range of 2-15 nm, while the primary carbon particles may have sizes of up to 80 nm based on the type of carbon support used. Early models of PEM fuel cells regarded CLs as interfaces with no consideration of their structural features. However, the significance of these structural features to the cell performance has been established. A particularly important issue that has resulted in significant efforts in CL modeling is the additional transport resistance observed at lowered Pt loading or with high loaded electrodes after degradation inducing cycles. 45,46 The experimental approach has utilized limiting current measurements with varying gas composition and/or pressure to separate the pressure dependent and pressure independent transport resistances. 47 The transport resistance in the CL is almost entirely independent of pressure and can be estimated with such limiting current measurements. Using this approach, it has been found that the CL transport resistance increases at lower Pt loadings [48][49][50][51][52][53] and this increase is strongly dependent on the available Pt area for reaction. Therefore, the resistance appears to be due to the increased flux near each active site at lower loadings. 1 Temperature sensitivities were used to determine the contributions of Knudsen diffusion and permeation through the ionomer thin film to the electrode transport resistance. 54 The ionomer thin film was found to be the dominant cause of transport resistance in the CL. More recently, the impact of carbon support and its porous structure on the local reactant and bulk protonic transport resistances have been highlighted. 2,55 Particularly, micro-pores with an opening smaller than 2 nm have been found to limit the reactant access to the Pt deposited inside the carbon pores. Despite such efforts, the root cause of the increased resistance remains largely unknown. 56 Several hypotheses have been made, but neither has been thoroughly validated.
Numerous models have been proposed to investigate the distribution of critical variables throughout the CLs and unveil the cause of increased transport resistance at lower loadings. The agglomerate model has been the most popular one for this purpose. In this model, the Pt particles are assumed to be dispersed on the primary carbon particles, many of which are assumed to aggregate during the fabrication process to form larger agglomerates covered by an ionomer thin film. The pore space in the CL is divided into two parts: the primary pores between carbon particles in each agglomerate, and the secondary pores formed between the agglomerates. Several variations of this model have been proposed where the intra-agglomerate space is either filled with water (i.e., water-filled agglomerates) 57,58 or ionomer (i.e., ionomer-filled agglomerates). 59,60 Initially, a wide range of agglomerate sizes (100-1000 nm) had been used and significant variations of the ionomer film thickness (10-100 nm) had been reported to match the experimental data. 61,62 Cetinbas et al. have developed a hybrid method for reconstruction of CL microstructure 63,64 and reported an agglomerate size distribution between 25 to 300 nm with most agglomerates having a radius in the range of 75 to 200 nm. 64 Furthermore, the upper limit of the modeling values for the agglomerate size and film thickness is not corroborated by microscopy studies. 65 Therefore, the validity of this structural picture has come under further scrutiny. In light of these experimental observations, some have argued that the agglomerates probably do not exist and have proposed homogeneous models for the electrode. 62 Others have continued to use the agglomerate models with agglomerate radii as small as 40 nm, 22,66 which is essentially the size of a carbon primary particle. Of particular interest is the work by Nissan researchers, 66 who showed that the conventional flooded-agglomerate model is not capable of reproducing experimental results with small agglomerate size and partial ionomer coverage. They modified the model to incorporate transport resistance near the electrochemical surface and showed that the modified model successfully predicted the experimental trends. Generally, more recent models rely on interfacial resistance at either the ionomer-gas or the Pt-ionomer interface or both to reproduce experimental transport resistance values. Jinnouchi et al. used molecular dynamics simulations to associate such resistance with a dense ionomer layer near the Pt surface. 67 Overall, attributing the additional resistance to interfacial phenomena has become increasingly common in the literature.
Despite its commonality, the interfacial resistance has not been experimentally verified. In fact, Liu et al. measured transport resistance in ionomer thin films and found no evidence of interfacial resistance when 3D diffusion was taken into account. 68 The uncertainty surrounding ionomer thin film properties, such as water uptake, 36,69-71 ionic conduction, 72 and gas permeation, 65 which can be significantly affected by confinement and substrate interactions, 73 has further contributed to the ambiguity of the source of this increased resistance. Some recent works have disputed the interfacial resistances or downplayed its significance. Darling has proposed an agglomerate model, in which the increased resistance is mostly attributed to the spherical diffusion through the agglomerate. 74 Others have investigated the inhomogeneity of mass fluxes near the Pt particles in agglomerates and the overlap between several agglomerates as possible culprits. 75,76 Most recently, Muzaffar et al. 77 have investigated literature data with a previously developed agglomerate model 78 and found that the reduction in Pt loading probably leads to higher levels of flooding in both the CL and GDL due to reduced vaporization capability of the CL with decreased thickness. They also elevated the fact that both experimental 79 and numerical studies 80,81 show only partial coverage of catalyst particles with ionomer, leaving an alternative transport path for oxygen to reach the active sites without facing the interfacial resistance at the Pt-ionomer interface. Therefore, they concluded that the increased transport resistance may be attributable to reduced oxygen diffusivity due to pore blocking effects of liquid water and the interfacial resistance does not play a significant role. The importance of water management in successful use of low loaded electrodes was also pointed out by Srouji et al. 82 The preceding literature review shows that the structural picture of the electrodes and the understanding of the factors that contribute to the transport resistance are still incomplete. Therefore, further model development and experimental investigations are required. Nevertheless, it should be noted that for the purpose of a full cell simulation, most of the proposed models can be parameterized to capture the local oxygen transport resistance, which is the most critical outcome of such models. Moreover, Kulikovsky has demonstrated that under certain conditions that are most relevant to typical fuel cell operation, the agglomerate model is not required. 83 Therefore, unless the goal of the model is to investigate different electrode designs at the nanoscale, a homogeneous model will be sufficient. Here, we use the model proposed by Hao et al., 62 which was shown to appropriately capture the increased resistance at lower loadings. The model achieves this by assuming full ionomer coverage and introducing rather significant interfacial resistances, which, in light of the above discussion, are disputable. Nevertheless, it is the general trend of the variations in the transport resistance that is required for our purposes. The model is briefly presented here and the reader is referred to the original publication for further details. 62 The model assumes Pt particles are deposited on primary carbon particles that are covered by an ionomer thin film. Liquid water in the pores of the electrode forms a thin film on top of the ionomer. This structural picture is used to derive the volume fraction of each phase (Pt, carbon, ionomer, and pore space) in both the anode and cathode CLs. However, the local transport resistance to hydrogen in the anode CL is assumed negligible and the calculations are only carried out for oxygen transport resistance. The oxygen in the pore space has to (1) dissolve in water, (2) diffuse through the water film, (3) dissolve in ionomer, (4) diffuse through the ionomer film, and (5) be adsorbed on the Pt surface. The model does not account for spherical diffusion, but uses instead a 1D diffusion equation to calculate the local flux of oxygen: [22] where N O 2 , c pore O 2 , and c Pt O 2 are the oxygen flux near the Pt surface, oxygen concentration in the CL pore space, and its concentration at the Pt surface, respectively. R T is the total local transport resistance: where the first, third, and last terms describe the interfacial resistances at the liquid film, ionomer film, and Pt surfaces, respectively. The fractional terms denote diffusional resistance through the water and ionomer thin films. A key argument made in developing the model is a geometrical one, where an effective diffusion length through the ionomer is calculated based on the effective surface area of a single Pt particle and the effective ionomer surface area available for that particle: [25] where r Pt and r c are the Pt and carbon primary particle radii, respectively, θ Pt denotes the fraction of Pt surface not covered with oxide species (see Reaction kinetics section), δ ion is the ionomer film thickness, and n Pt is the number of Pt particles deposited on a single carbon particle. The effective ionomer film thickness is then calculated by: The same scaling factor is used to scale the interfacial resistance at the Pt surface: This scaling is one of the most important features of the model, as it compensates for the fact that 3D spherical diffusion is neglected, and allows for the effects of high fluxes near sparsely deposited Pt particles to be captured by the model. It is imperative, however, to be cautious and not put too much emphasis on the source of the local transport resistance in this model. As mentioned earlier, the electrode structure assumed in this model is contentious. Nevertheless, on a macro-level, the predictions match the experimental observations, which is the most important aspect for full cell simulations. Finally, another important assumption made in the model is that the interfacial resistances are proportional to the diffusional resistances. This is done due the lack of measured data for the interfacial resistances at various interfaces. In particular, three fitting parameters k 1 , k 2 , and k 3 are introduced: Therefore, the various terms contributing to the transport resistance are identified. Noting that: [29] where i ORR is the volumetric ORR current density, a c is the volumetric surface area of the ionomer, and x is the number fraction of carbon supported Pt particles (used to model the effects of catalyst dilution by bare carbon), Equation 22 can be written as: This algebraic equation can be solved numerically to find the oxygen concentration at the Pt surface. It is worth pointing out that an analytical solution is possible in the case that reaction order is assumed to be unity for ORR. 62 Nevertheless, such an assumption may be unrealistic and in some cases inconsistent with the ORR kinetics model (see Reaction kinetics section). Therefore, we use the numerical solution with no assumption on the reaction order for ORR to avoid such inconsistencies. It is also important to have a consistent set of structural parameters for the CLs. In particular, volume fraction of different phases ought to be known. These volume fractions can be calculated as follows: 60,63 where ε i is the volume fraction of i, L c/PT is the carbon/Pt loading, ρ i is the density of i, δ CL is the CL thickness, and (I/C) denotes the ionomer to carbon ratio. The remaining CL volume constitutes its pore space (ε CL = 1 − ε c − ε Pt − ε ion ). Finally, the ionomer and liquid water film thicknesses are given by: . [35] This completes the CL model used in this work. The reader is referred to 62 for further details about this model. As for the parameter values, an I/C ratio of 1.1 and an electrochemically active area (ECSA) of 65 m 2 Pt /gr Pt are assumed for both anode and cathode CLs. The anode Pt loading is assumed to be 0.1 mg/cm 2 with a Pt/C weight percentage of 30%, while the cathode Pt loading is changed between 0.4 and 0.05 mg/cm 2 , considering a Pt/C weight percentage of 40% in all cases.
Reaction kinetics.-Accurate models for the HOR and ORR half reactions are required for the model. HOR is known to have facile kinetics and does not result in significant performance loss under most typical conditions. Therefore, it is typically described using a simplified Butler-Volmer kinetics model. Here, we use the dual-pathway kinetics model proposed by Wang et al., 84 where the volumetric current density is found by: where a Pt is the active volumetric surface area of Pt, i 0T and i 0H are the exchange current densities for the Tafel and Heyrovsky pathways, respectively, ϑ is a potential constant, and η HOR is the anode overpotential.
The ORR kinetics are more complicated than the HOR and require further attention. Again, various forms of the Butler-Volmer model have been used to describe the ORR kinetics. More recently, the effects of surface coverage have been considered to derive a modified Tafel expression: 85 where i 0,ca , α ca , η ORR , γ ca , are the cathode exchange current density, transfer coefficient, ORR overpotential and reaction order, respectively, and ω denotes the energy parameter for the Temkin isotherm. The model results in a potential dependent Tafel slope. The oxide coverage is potential and time dependent as cyclic voltammograms (CV) show considerable difference between the anodic and cathodic sweeps. 85 A simple sigmoidal curve is fitted to steady-state measurements: 62 , [38] where E is the cathode potential vs. reference hydrogen electrode (RHE). A more elaborate model for ORR kinetics is the double trap (DT) model originally proposed by Wang et al. 86,87 The model includes two pathways for oxygen adsorption: a reductive adsorption (RA) and a dissociative adsorption (DA) pathway. The latter is followed by a reductive transition (RT) to adsorbed OH. In either case, the adsorbed OH is desorbed through a reductive step (RD) to form water. The original formulation neglected the reverse RD step and concluded that ORR activity is limited by the desorption of strongly adsorbed O and OH. Moore et al. 88 modified the model by including the backward reactions and refitting the parameters and found ORR to be adsorption limited. Moreover, the coverage of adsorbed species predicted by the modified model tends to zero at high overpotentials, whereas a constant nonzero value was predicted with the original model. 86 The modified model is in better agreement with the experimental coverage values reported by Subramanian et al. 85 Other modifications to the DT model have been proposed as well. Markiewicz et al. 89 added two elementary reactions to the model: a reductive addition of a proton to oxygen molecule, producing an adsorbed protonated superoxide, and another reductive addition of proton followed by dissociation into adsorbed OH. Through these modifications, they reported a significant coverage of Pt sites by adsorbed HO 2 species at high overpotentials. More recently, Jayasankar et al. 90 replaced the DA step with an associative adsorption (AA) into adsorbed HO 2 , which is followed by dissociative transition steps into adsorbed O and OH. They have also extended the model to include oxide growth mechanisms. Their results corroborate those of Markiewicz et al., as they also find an increase in HO 2 coverage at high overpotentials. This can have significance for studies with low loaded catalysts, as it provides another possible explanation for the reduced performance observed experimentally.
In this work, we use the modified DT model proposed by Moore, 88 as it has been parameterized for fuel cell polarization curves and used by others in full cell simulation. 19,22 In this model, the ORR current can be described as the current from a single RD step: where i * is a reference prefactor (similar to the exchange current density in the Butler-Volmer model), k is the Boltzmann constant, G * RD and G * −RD are the potential dependent activation energies of the forward and backward RD step, respectively, and θ i denotes the coverage of species i. The expressions for the activation energies and species coverage can be found in Ref. 88.
The DT model is used for simulation case studies. However, when comparing with experimental data, we have chosen to work with the Tafel model in Equation 37, as its parameters are more intuitive and allow for easier parameterization of the model and can also reproduce the kinetic current predicted by the DT model with a varying reaction order. 91 It should be noted that several effects have been neglected to simplify the model and avoid ambiguity in the results. First, the steadystate coverage profiles are used in the kinetic equations and the dynamics of oxide growth are ignored. These dynamics can be very slow as observed in low frequency impedance spectra 92,93 and coulometric measurements. 94 Such dynamics can result in a hysteresis loop in the Tafel plot obtained through CVs even when a low potential prehold is used to reduce the oxide layer. 95 Therefore, oxide growth dynamics can have a profound impact on current transients, especially at higher potentials. However, including these dynamics adds to the complexity of a model, whose main focus is on mass transport and hydration effects. Hence, the oxide growth dynamics are neglected in this work. It should also be pointed out that the ORR activity is shown to be affected by presence of ionomer. [96][97][98] This effect is not explicitly taken into account in the current model, since doing so will add to the uncertainty in the parameter set. Nevertheless, the exchange current density values (or the reference prefactor in the case of the DT model) used are supposed to capture this reduced activity.
Finally, the dependence of ORR kinetics on the relative humidity (RH) is also neglected in this work. This effect was reported by Xu et al. 99 to be significant, resulting in up to 100 mV difference under dry conditions, even when protonic resistance in the CL was taken into account. 100 However, work by GM shows much less pronounced effects of RH on ORR kinetics. 101,102 This discrepancy in the reported values could also be partly due to the effects of RH on water oxidization and subsequent catalyst poisoning. 103 The accessibility of Pt in the inner pores of porous carbon support is also shown to decrease at low RH values, which can result in loss of electrochemically active area. 104 Regardless, the RH effects on ORR kinetics may be included in the model by scaling the exchange current density in the Butler-Volmer model (i 0 ) or the reference prefactor in the DT model (i * ). A scaling factor varying linearly with ionomer water content has been used for this purpose by Gerteisen et al. 9 More recently, a scaling factor that changes with ionomer water content in a sigmoidal fashion has been proposed, 42 which is in a better agreement with the experimental trends. Such scaling factors may be treated as fitting parameters in performance models to enhance the predictive capabilities. Nevertheless, we have chosen to leave this factor out in order to simplify the model and allow for a clearer understanding of the transport phenomena.

Mixed wettability model for porous layers.-
The main goal of a model for the porous layers is to define a mapping from operating conditions and material properties to effective charge, heat, and mass transport properties. This problem has been studied on a variety of length scales ranging from microscopic lattice Boltzmann 16 and pore network modeling studies 105 to macroscopic models with empirical relationships. 106 The microscopic models, along with significant advances in experimental techniques to characterize porous layers at higher resolutions, can be used to develop a fundamental understanding of various transport phenomena in such layers. Even though such models cannot be used in full cell simulations due to significant computational requirements, they can be utilized to refine the macroscopic models of lower complexity.
Understanding the water phase change process and its transport through the porous layers is also of crucial importance. To this end, one particular model for porous layers that has gained more popularity in recent years is the mixed wettability pore size distribution (PSD) model that was proposed by Weber et al. 107 The model represents the pores as bundles of cylindrical capillaries that are randomly joined together using log-normal distributions. The key feature of the model is that it accounts for mixed wettability of the layers, which is ignored for the most part in many of the macroscopic models. Therefore, both hydrophilic (HI) and hydrophobic (HO) pores are considered to derive PSDs and contact angles. The original implementation by Weber et al. 107 assumed the HI and HO PSDs to be identical. Furthermore, a two-point discrete contact angle distribution was assumed. A similar model was used by Eikerling for transport studies in the cathode CL, although he did not consider mixed wettability, choosing to investigate the PSDs due to primary and secondary pores in the CL. 108 More recently, Villanueva studied effects of different PSDs for the HI and HO pores. 20 However, recent implementation of the model in a full cell simulation by the same group seems to be using similar PSDs for both HI and HO pores. 19 It is worth mentioning that this model was further extended by Weber to include a continuous contact angle distribution (CAD). 109 This extension was shown to improve the predictive capabilities as well as the numerical robustness of the model for use in full cell simulations. A continuous CAD with a discrete PSD was used by Cheung et al. 110 Nevertheless, adoption of the continuous CAD has remained minimal in the literature due to unavailability of CAD for most of the porous layers of interest.
This work utilizes the mixed wettability model with an implementation that allows for different PSDs to be used for HI and HO pores and also includes the continuous CAD. However, we obtain PSDs and contact angles from the literature, and therefore our implementation coincides with the original implementation by Weber et al. 107 when continuous CAD is not available.
The model equations can be found in the literature 20,109 and are omitted here for space considerations. The inputs to the model include the PSDs, the fraction of HI pores, and CADs (currently two-point discrete CADs are used). The model is used to obtain water retention curves (liquid saturation vs. capillary pressure), relative permeabilities of the gas and liquid phases, Knudsen radii, and liquid-gas interfacial area available for phase change in the CLs, MPLs, and GDLs. The model calculations are conducted off-line and the resulting curves are used in the full cell simulations using cubic-spline fitting. The model parameters used for the simulation case studies in this work are presented in Table IV and the resulting water retention curves are shown in Fig. 2. The CL parameters used in this study are those reported by Mashio et al., 111 who obtained experimental PSDs through nitrogen adsorption. In particular, the PSD for a CL with graphitized Ketjen Black carbon support and an ionomer to carbon ratio of 0.9  19,23 for SGL 34 series that are also applicable for the SGL 24 series used in this work.
It should be noted that, as has been shown by Zenyuk et al., 17 the GDL PSD changes with compression. Therefore, it seems reasonable to use a PSD corresponding to higher compression under the land compared to the one used under the channel. However, using the PSDs that were reported for SGL series by Zenyuk et al. under various compressions, 17 we found that the effective transport properties of interest show little change with the PSD variations at the compression levels of interest (1 to 1.5 MPa). Therefore, the changes in the PSD with compression under the land are ignored in this work. Finally, as it has been alluded to by Weber, 109 the bundle of capillaries model breaks down for wide PSDs, which in turn results in very low relative permeabilities predicted by the model. We have found this to be especially problematic for the MPL. Therefore, a 5-th order power law is used to estimate the liquid and gas phase relative permeabilities for the MPL.
Calculation of effective properties.-To complete the model formulation, effective properties such as gas diffusivity and thermal conductivity values are needed. Some of the layers demonstrate rather considerable anisotropy due to their heterogeneous structure, which should be taken into account. Furthermore, the effects of nonuniform compression under the channel and lands should be considered to obtain an accurate in-plane distribution of the variables of interest. 40,112 Accordingly, we have carefully examined the literature for the reported values of such transport properties. When applicable, the land-channel variations in parameter values are applied in a continuous fashion using sigmoid functions. This is in better agreement with the observed pressure distribution and also simplifies numerical convergence.
In this work, we use SGL 24BC and Nafion 211 as the diffusion media and membrane, respectively. These materials are chosen due to their standard application in the fuel cell literature and an abundance of experimental characterization data available for them. The layer thickness and porosities are listed in Table V. Note that a compressed GDL thickness is assumed based on a compressive load of 1 MPa, which is expected to result in a strain of about 0.2. 115 While a uniform thickness is used for both the channel and land locations, the collapse of pore space is applied to the land area, where a reduced porosity of 0.69 is used for the GDL. The CL and MPL are assumed to be incompressible. Furthermore, note that Nafion 211 has no reinforcement, yielding ε ion = 1 in the membrane region. Finally, it should be pointed out that an intermediate composite region is believed to exist between the MPL and GDL with transport properties that are considerably different from those of either layers. Since the properties of this intermediate region are not well known, it is not explicitly modeled in this work. The material property variations between adjacent layers are taken into account using smooth sigmoid functions to improve convergence. A detailed discussion of the effective transport properties used in the simulation studies follows.
Effective diffusivity.-In calculating the diffusivity of species i, contributions from both molecular and Knudsen diffusion are taken into account: where D Kn,i is the Knudsen diffusivity and D mix,i is the molecular diffusion coefficient. Knudsen diffusivity is given by: where r Kn,i is the Knudsen radius of the porous layer, which is obtained from the mixed wettability model in this work, and M i is the molecular mass of species i. The molecular diffusion coefficient is given by: 116 where x j is the molar fraction of species j and D i, j denotes the binary diffusion coefficient of species i in j. 117 With D i available, the effective diffusivity is calculated as which accounts for the tortuous pathway for gas transport inside the porous layers as well as the pore blocking effects of liquid water accumulation. Several microstructure-property functional relationships have been proposed for both f (ε) and g(s) in the literature, most of which take the form of a power law. 118 f (ε) IP = 1.074ε − 0.335, [44] f (ε) TP = 0.906ε − 0.252, [45] where the subscripts IP and TP stand for the in-plane and throughplane directions, respectively. These relationships were suggested for SGL 25BA series, which do not have the MPL coating. Due to lack of data, the same relationships are used for the MPL. The reader should be cautious in applying these relationships to other types of diffusion layers, as they are explicitly derived for SGL carbon papers. Typical power laws are better applicable in general and are suggested for different types of diffusion layers. As for g(s), the following relationships are used: 118 g(s) TP = (1 − s) 2.15 . [47] These relationships were determined for Toray carbon papers and are used here due to lack of data for SGL series. It should be noted that the nearly isotropic relationships were developed for local conditions, i.e., domains with flat saturation distributions. 122 In contrast, Niu et al. found a more significant difference between the liquid saturation effects on the in-plane and through-plane diffusion coefficients, fitting the results with cubic and quadratic power laws, respectively. 119 Therefore, such functional relationships should not be taken for granted. Rather, we believe that it is a better practice to leave the order of dependence as a fitting parameter when experimental performance data are available. Finally, the correction factor for effective diffusivity calculations in the CL is calculated as follows: 19,108  where ε p is the percolation threshold, which is assumed to be 0.25 in this work, and H is the Heaviside function.
Absolute permeability.-A range of values for absolute gas and liquid permeabilities are reported in the literature. In most of the cases, the absolute permeability of the GDL is found to be on the order of 10 −11 m 2 . 123 The MPL permeability values are typically one to two orders of magnitude smaller than those for the GDL. 124 In spite of such measurements, it has been shown that these permeability values will result in a negligible pressure drop across the porous layers due to oversimplification of the capillary dominated transport through the use of Darcy's law. 22,125 Therefore, it has been suggested that the experimentally reported values should be reduced by several orders of magnitude to obtain a realistic pressure drop. 22 In addition to these considerations, Holzer et al. recently found that the through-plane permeability values are slightly higher than the in-plane values for SGL 25BA GDLs. 123 Taking these into account, the absolute permeabilities assumed in the model are given in Table VI.
Thermal conductivity and heat capacity.-There is an extensive literature on the thermal conductivity of the PEM fuel cell layers through both modeling and experimental means with somewhat scattered results. In selecting the thermal transport parameters, one has to pay attention to the changes in thermal conductivity with liquid saturation, in addition to the anisotropy and compression effects. Another difficulty is in distinguishing between the thermal properties of the MPL and the GDL from the data obtained from a composite layer. Here we briefly review the existing literature for SGL papers.
One of the earliest works in this area was by Khandelwal et al., who measured the TP thermal conductivities of various cell layers. 126 They reported a value of 0.31 W/(m · K) for SGL BA series (with no MPL). Unfortunately, they did not report the number specification of the GDL. This can bear some significance as the SGL 24 series have more binder that can improve the fiber to fiber contact and increase the thermal conductivity. 127 Nevertheless, this value is well within the range of 0.26-0.37 W/(m · K) reported by others for the same type of GDL. [127][128][129][130] Accordingly, the base value of GDL TP thermal conductivity is set to 0.3 W/(m · K) under the channel and to 0.45 W/(m · K) under the land due to the inhomogeneous compression. 130 As for the IP thermal conductivity, a base value of 12 W/(m · K) is used 131,132 for both the channel and land locations as the effect of compression on IP conductivity is assumed to be minimal.
For the MPL thermal conductivity, the reported values are far more inconsistent than those for the GDL. Such discrepancies stem mostly from unknown contact resistances, uncertainties about the MPL thickness in a combined layer, assumed compressibility or incompressibility of the MPL with applied pressure, and the nature of the transition region between the MPL and GDL. These have resulted in reported values for the TP thermal conductivity ranging from 0.035 133 to 0.6 W/(m · K). 130 An interesting observation was made by Burheim et al., 134 who argued that the MPL has a lower thermal conductivity than the GDL (0.08 W/(m · K)), with an intermediate composite region between the two layers that has the highest thermal conductivity with an essentially flat temperature distribution. Here we use a value of 0.15 W/(m · K) for both the channel and land locations. This is based on the assumption of incompressibility of the MPL, which has been questioned recently. 134,135 Nevertheless, this value is in the range of reported values in the literature. The base value for the IP thermal conductivity of the MPL is chosen to be 3 W/(m · K) based on the literature. 131 The reported TP thermal conductivities for the CL range from 0.04 136 to 0.34 W/(m · K). 137 In this work, we use the base value of 0.27 W/(m · K) reported by Khandelwal et al. 126 for both the IP and TP thermal conductivities assuming no anisotropy for the CL. 137 Liquid accumulation in the pores can alter the thermal conductivity of the porous layers. In this work, we use the following approximation to capture this effect for the TP thermal conductivity of the GDL: 115 k eff T = k T,base + 1.44s, [49] where k T,base is the base value reported in Table VI. For the TP thermal conductivity of other layers (CL and MPL) and the IP conductivity of all porous layers, volume averaging is employed: where k T,l = 0.569W/(m · K) is the thermal conductivity of liquid water. Finally, the thermal conductivity of the membrane in both the IP and TP directions is given by: 138 The volumetric specific heat capacities (ρc p ) used in the model are: 1.9, 1.562, 10  [52] The debated suppression of ionic conductivity in thin ionomer films is not taken into account in this work. 142 Future work should aim at investigating such effects through parametric studies.
Evaporation and condensation rates.-A value of 2 × 10 −2 mol/(cm 2 · s) is used for the condensation rate to avoid the nonphysical case of oversaturated gas phase. 19 The evaporation rate is set to 2 × 10 −3 mol/(cm 2 · s). Even though this value is lower than the condensation rate, it yields rather fast evaporation kinetics, which agrees with the experimental findings of Zenyuk et al. 143 that showed the evaporation to be transport-limited. The discrepancy between the evaporation and condensation rates is corroborated by experimental findings in the literature. 144 Furthermore, the rate of phase change is  144 which is not taken into account in this work. It should be noted that the evaporation rate is a critical parameter and may have a significant impact on water balance in the cell depending on the operating conditions. Future work should aim at a sensitivity analysis for this parameter to better understand its impact on the performance in the two phase regime.
Boundary conditions.-The model boundary conditions (BCs) are given in Table VII, where n denotes the unit normal vector. Symmetry boundary conditions (i.e., zero flux) are applied at the top and bottom boundaries of the modeling domain shown in Fig. 1. The temperature BCs include two heat flux BCs at the channel and land locations. The channel heat flux corresponds to convective heat transport with the gas stream (h conv = 0.2 W cm 2 ·K ), while the land BC accounts for the thermal contact resistance (R T,cont = 2 cm 2 ·K W ) between the plate and the GDL. The molar fractions of gas species are also modeled with mass flux BCs at the channel location to account for the convective mass transport resistance. In the corresponding equation, D free i, j denotes the bulk diffusivity of species i in species j, Sh is the dimensionless Sherwood number (=2.7), and D h is the hydraulic diameter of the channel. Dirichlet BCs are used for gas pressures at the channel boundaries. The liquid pressure BC requires further attention. Various types of BCs have been used for this purpose, including Dirichlet BC for liquid saturation or capillary pressure, 145 as well as Neumann type BC. 9 In this work, we use the following BC: where k l,flux is a parameter determining rate of water outflow, s 0 controls the liquid saturation at which water outflow begins, and σ s is a dimensionless parameter used to smooth the transition between no flux BC and the outflow BC. Note that the parameter s 0 essentially accounts for the break-through pressure, which is the capillary pressure required for liquid water to flow out of the porous GDL. The values of the three parameters used in this work are: k l,flux = 8 × 10 −4 g cm 2 ·s , s 0 = 0.1, and σ s = 0.01. It should be pointed out that this BC can be parameterized to be identical to the BC used by Zhou et al. 19 However, it has the advantage that the parameters are more intuitive, which can simplify the parameterization process.

Numerical implementation and model validation.-The model
is implemented in the commercial finite element software COMSOL Multiphysics 5.3a. A mapped mesh consisting of 5080 quadrilateral elements is used throughout the domain with increased mesh density in the membrane and catalyst layers. Furthermore, the mesh density is exponentially increased near the boundaries between adjacent layers to accommodate the different material properties. The backward differentiation formula (BDF) method is used for time stepping and the maximum time step size is limited to 200 milliseconds. The resulting linear system is solved using the MUMPS direct solver provided in COM-SOL. To improve the computational efficiency, an under-relaxation scheme is employed, where the value of liquid saturation at the pre-vious time step is used to calculate effective properties such as the diffusion coefficients at the current time step. This was achieved using the Previous Solution operator in COMSOL 5.3a. In a preliminary study, it was found that the under-relaxation scheme can result in up to five times faster solutions in the two-phase regime. The results for the 340 seconds long simulation case studies in this paper were computed in 5 to 12 hours depending on the condition, with the most difficult cases being the ones where the transition from dry to wet conditions takes a long time. The simulations were run on a desktop computer with a 3.5 GHz processor and 16 GB of RAM.
The model is validated with experimental data by Gerteisen et al. 9 The results can be found in the Supplementary Information accompanying this paper.

Simulation Case Studies
To better understand the transient behavior of the cell, several simulations are conducted using the model developed in this work. In particular, the transient performance under a variety of temperature and humidity conditions as well as different Pt loadings in the cathode CL is investigated. Furthermore, we investigate the cell dynamics under both potential and current control operating modes. The former constitutes running the model with voltage as an input, while the latter takes the cell current density as the input. As shown below, the dynamics of the cell response can be dramatically different depending on the operating condition. All of the simulations in this work were conducted at a pressure of 1.5 bar for both sides. The gas feeds are assumed to be pure hydrogen and air for the anode and cathode sides, respectively. Finally, same RH values are used for both the anode and cathode sides and initial conditions for all simulations are identical.
Potentio-dynamic simulations.-The first set of simulations are those under voltage control or potentio-dynamic mode of operation. Here a voltage profile is applied and the current density is allowed to vary with time. The time varying current density also means that the rate of water production changes with time, which complicates the analysis of the dynamics to some extent. Nevertheless, useful insights can be obtained from these simulations.
The voltage profile for these simulations is shown in Fig. 3. The profile is made up of the following voltage steps: 0.8-0.6 V, 0.6-0.4 V, 0.4-0.6 V, and 0.6-0.8 V. Note that the step changes are smooth and happen over a period of 1 second for numerical convergence. This profile allows us to inspect the transients during both load increments and decrements. The 100 second hold time used at 0.6 and 0.4 V does not allow the system to fully reach its steady state conditions. Nevertheless, this hold time is limited due to computational reasons and is long enough for the model to settle to a quasi steady state before another change in the load.
Overall, 36 simulations are conducted under the potentio-dynamic mode based on a full factorial design for variations in RH (30, 60, and 90%), operating temperature (40, 60, and 80°C), and cathode Pt loading (0.4, 0.2, 0.1, and 0.05 mg/cm 2 ). Note that the CL thickness is assumed to scale linearly with the Pt loading. The resulting current   Based on these results, several conclusions can be made about the through-the-membrane phenomena affecting the transient response of the cell. The following analysis of the average response is organized based on the step change in the load. Discussions on the distribution of the critical variables are provided later in the paper.
Voltage step from 0.8 to 0.6 V.-During this step change, the current responds monotonically with varying settling times that decrease with channel RH, i.e., a faster response for more humidified conditions, which can be attributed to sufficient membrane humidification. This can be observed in Fig. 5, which shows the membrane water content dynamics, where the transient response is found to significantly depend on the operating conditions. In particular, under dry conditions that dry out the membrane prior to the voltage step down, the membrane water content increases monotonically with the step change in voltage. This increase is less pronounced at higher temperatures, where further increase in temperature at higher loads results in lower water uptake by the membrane. As the humidity increases and the membrane holds enough water in its initial state prior to the step change, we observe some cases with reverse response, i.e., an initial decrease in the membrane water content followed by an increasing trend (see, for example, the case with T = 40°C, RH = 60% and a Pt loading of 0.4 mg/cm 2 in Fig. 5). This reverse response is due to EOD that tends to dry out the anode side of the membrane and is only seen at higher current densities. Another observation is that the changes in the membrane hydration are much more pronounced at lower temperature, where slight variations in water production rate can significantly alter the membrane water content. The slow relaxation dynamics discussed in the Model Formulation section are also evident in Fig. 5.  In particular, we note that the relaxation dynamics become slower at higher water contents. These relaxation effects are not observable in the current dynamics, since the ohmic drop at this relatively low load is insignificant. Finally, as shown in Fig. 6, some liquid water builds up in the cathode GDL after the voltage step down for three conditions with low temperature and high humidity. As expected, the build up of liquid is faster under cooler and wetter conditions, while for some conditions the dynamics are slow and the liquid saturation does not reach a steady state within the 100 seconds hold at 0.6 V (see, for example, the case with T = 40°C, RH = 60% and a Pt loading of 0.4 mg/cm 2 in Fig. 6).
The figures also demonstrate that the cathode Pt loading has an impact on the transient response by influencing the current generation, membrane hydration, and liquid saturation dynamics. More specifically, Fig. 4 shows that higher current densities are achieved with higher loadings, which in turn affect the membrane humidification process, especially under drier conditions, where the membrane is humidified with the electrochemically generated water. This can be clearly seen in Fig. 5 for T = 40°C and RH = 30%. It is observed that with higher loadings, the water generation is high enough to humidify the membrane, whereas with a loading of 0.05 mg/cm 2 , the membrane remains dry after the step change in the voltage. Additionally, Fig. 5 shows that the membrane water content is indeed influenced by the Pt loading, which can be associated with different levels of current and heat generation as well as a change in the overall water balance in the cell with the CL thickness. It can also be observed that the cases with higher cathode Pt loading show higher levels of liquid saturation in the GDL (Fig. 6). This can be attributed to the fact that higher Pt loading results in higher current density and therefore higher rate of water generation. In addition, the resulting variations in the CL thickness with changes in Pt loading mean that lower loaded CLs generate more heat on a volumetric basis, which creates a stronger drive for water evaporation. The impact of Pt loading on the overall water balance in the cell is discussed in further detail below as well as in the discussion of the galvano-dynamic simulations, where the results are not convoluted by varying levels of water generation.
Voltage step from 0.6 to 0.4 V.-The second voltage step results in more involved dynamics in some cases, but the current response can be categorically identified as being monotonically increasing or displaying an overshoot. In particular, the drier conditions tend to result in a monotonic increase in the current density (see the first row of Fig.  4). Similar to the previous step, this monotonic response is associated with a hydrating membrane. This is mostly evident, for instance, at T = 60°C and RH = 30%, where the relaxation dynamics for the membrane water uptake also play a role in the slow increase in current density. Under wetter conditions, however, the gas phase in the cathode CL is saturated with vapor and the membrane protonic resistance is low enough prior to the step change. This low protonic resistance can support high current generation immediately after the step change. The high current density dries out the anode side of the membrane with EOD and increases the protonic resistance, which results in a performance drop as seen in Fig. 4. The overshoot response due to EOD is relatively fast and settles within 5 seconds of the step change, when the generated water on the cathode side diffuses back toward the anode and rehydrates the dry portion of the membrane. 146 It is also observed that the overshoot becomes progressively less significant as the temperature increases. The large overshoots at low temperatures can be attributed to the high sensitivity of the membrane hydration state to changes in the current density. This high sensitivity stems from more rapid changes in the environmental conditions (T and RH) in the CL and the increased protonic resistance at these lower temperatures. It should be noted that flooding of the porous layers, discharging of the electrochemical double layer, and mass transport limitations are also believed to lead to this type of behavior. 4 However, the overshoot in the presented results is due to membrane dry out with EOD, which is in agreement with other experimental results. 146 Following this overshoot, the wetter conditions display a relatively slow drop in performance during the hold at 0.4 V (e.g., see the case with T = 60°C and RH = 90% in Fig. 4). This slow drop in performance is attributable to liquid build up in the GDL that incurs a mass transport resistance and is most significant for the cases where the current density goes above 1.5 A/cm 2 .
The case with T = 60°C, RH = 60%, and a Pt loading of 0.4 mg/cm 2 displays interesting dynamics after the voltage step. An initial overshoot due to EOD is observed that causes the current density to drop by about 0.1 A/cm 2 within 2 seconds of the step change. This drop is then followed by an increase of about 0.06 A/cm 2 during the next 15 seconds. This increase is due to liquid accumulation in the CL that helps the membrane humidification without causing mass transport limitations. Afterwards, the slow current decay due to liquid accumulation in the cathode GDL can be observed, which continues until the next voltage step. In fact, this particular order of liquid build up in the porous layers (the CL pores followed by those of the MPL and GDL) is seen under most typical conditions. However, the observability of this behavior from measurements of current alone depends on the water retention capabilities of the different layers as well as the operating conditions used for the experiments.
As for the membrane water content, the most notable observation is that for drier conditions where the gas phase in the CL remains unsaturated, this voltage step results in better membrane humidification due to higher rates of water generation (Fig. 5). However, this trend is reversed for wetter conditions, where the water content drops after this voltage step. This drop is again attributable to water removal to the cathode through higher EOD at higher current densities. Additionally, the variations in membrane water content with the voltage step are more significant at lower temperatures as was the case during the previous step change. Average liquid saturation in the cathode GDL also exhibits trends similar to those for the previous step (Fig. 6).
Voltage step from 0.4 to 0.6 V and from 0.6 to 0.8 V.-Similar to the previous steps, the current response to the step increase in voltage is either monotonically decreasing (as is the case at T = 60°C and RH = 30%) or exhibits an undershoot. Again, this behavior can be directly correlated with membrane water content and further discussion is omitted here. Instead, we focus on the trends during the dry-out phase when the voltage is increased. In particular, we note that under some drier conditions, the performance starts to decay after a while. This is seen, for instance, for all Pt loadings at T = 40°C and RH = 30% in Fig. 4. Looking at the membrane water content dynamics in Fig. 5, we note that this decay in performance is directly related to membrane water loss. This behavior can be explained by a moving evaporation front that starts in the GDL and progresses toward the CL as time goes on. Therefore, immediately after the step change the ionomer in the cathode CL is in contact with a liquid reservoir, which improves water uptake. As the evaporation front reaches the CL and the accumulated liquid evaporates, the membrane starts to lose water, which results in further performance decay. In these simulations, the time delay between the step change in voltage and the evaporation of CL liquid water depends on the operating conditions as well as the CL thickness, with hotter conditions and thinner CLs generally resulting in the shortest time delays. More generally though, this time delay is determined by the HI contact angle of the CL as well as the evaporation rate used in the model. A lower HI contact angle makes evaporation of water in the HI pores more difficult and prolongs the time delay, whereas a high evaporation rate reduces this delay.

F3167
Finally, it should be mentioned that based on experimental results, a better performance may be expected at 0.6 V after the hold at 0.4 V (220-320 seconds) compared to the performance at 0.6 V before the hold at 0.4 V (20-120 seconds). This performance gain is attributable to better membrane hydration as well as clearing of the Pt sites from oxide coverage at low potentials. This improved performance will diminish slowly as the membrane dehydrates and oxide species grow on the Pt surface again. As explained in the Model Formulation section though, this work only captures the former dynamics, as the oxide growth is ignored in our model and steady state coverage values are used.
The characteristic responses observed after the last voltage step are similar to those discussed so far and further discussion is omitted here.
To further investigate the transient phenomena through the cell's thickness, Fig. 7 illustrates the average liquid saturation in the CL and GDL under both channel (CH) and land (LN) regions of the cell, the average ionomer water content in the anode and cathode CLs, and the normalized membrane water flux defined as: , [54] which is averaged over the area of the membrane. As a convention, a positive value denotes water flux toward the cathode. The presented results are for the cold and dry conditions (T = 40°C and RH = 30 %) with high (0.4 mg/cm 2 ) and low (0.05 mg/cm 2 ) Pt loading in the cathode CL. First, we observe that immediately after the second step decrease in voltage (0.6 to 0.4 V at 120 seconds), the cathode CL becomes flooded. This flooding takes place within 10 seconds of the step change. Also, note that the CL flooding occurs slightly faster under the land location compared to the channel location. After all the hydrophilic pores in the CL are filled, water starts to condense in the GDL. Most of the condensation happens under the land, where lower temperature and higher resistance to vapor transport promotes the phase change process. This liquid water then flows toward the channel location. This can be seen in the figure, as there is a delay between liquid accumulation under the land and channel regions of the GDL. This delay is governed by the time it takes for the liquid water condensed under the land to reach the channel location. After the voltage is increased back to 0.6 V at 220 seconds, we see that the GDL dry out is initiated under the channel. The dry out happens at a slower pace under the land location. Once the GDL is completely dry, the CL starts to lose its liquid water.
As for the ionomer water content in the CL, two main observations can be made. First, at lower loads the ionomer water contents in both CLs are close and as the load is increased, a more significant distribution develops across the membrane thickness with intensified EOD and back diffusion. Second, we note that after the load is decreased, the CLs maintain a high ionomer water content as long as the liquid water in the CL has not evaporated. After that liquid water has evaporated though, the cathode ionomer loses water to its pore space, which also diminishes the water back diffusion to anode, which in turn results in dry out of the anode CL. This behavior can also be seen in the last column of Fig. 7, where the smallest values for β are obtained when the cathode CL has liquid water while the anode CL is dry, which is in agreement with experimental measurements by Adachi et al. 147 It is also seen that this flux is significantly reduced after the liquid reservoir in the cathode CL has evaporated (see the plots between 220 and 300 seconds). This is specially pronounced for the case with low Pt loading.
Another important observation from the normalized water flux plots is the significant overshoots and undershoots during the step changes. This behavior is associated with EOD that immediately drives water to the cathode side, whereas the back diffusion requires time to establish a balancing water flux to counter EOD. Such transients qualitatively agree with experimental measurements. 148,149 In addition, it can be seen that higher current densities generally tend to force more water toward the cathode (larger β values) due to intensified EOD, which is in agreement with experimental results. 150,151 Finally, the distribution of critical variables for some conditions are shown in Figs. 8-9. In particular, the distributions for the membrane and cathode temperatures, cathode liquid saturation, ionomer water content, and volumetric ORR current density are shown for high (0.4 mg/cm 2 ) and low (0.05 mg/cm 2 ) Pt loadings. These distributions are obtained at the end of the hold at 0.4 V. The temperature plots show a rather significant temperature gradient (about 1.5°C) across the MPL thickness, which is due to its low thermal conductivity. This low conductivity also results in heating the CL and enhances water evaporation. 23 The cathode MPL and GDL remain free of liquid water under these hot conditions, whereas the hydrophilic pores in the cathode CL are filled with liquid water for the hot and wet (T = 80°C, RH = 90%) conditions. A considerable gradient of water content is established across the thickness of the CCM, with a dry anode CL and a wet cathode CL. There is a close correspondence between the location of maximum water content in the cathode CL and the volumetric rate of ORR. In particular, under the dry conditions (Fig. 8) resistance is a major contributor to performance loss. Therefore, the highest volumetric current is observed under the land location, where the membrane water content is highest. At higher humidities (Fig. 9) the location of maximum current generation moves toward the channel region, where the mass transport limitations are minimal. Furthermore, a higher portion of the Pt is utilized, as the region close to the MPL is not severely limited by proton transport resistance. Comparing the low and high Pt loading cases, the most notable difference is in the volumetric current distributions, which stems from the thinness of the CL with low Pt loading and the resulting increase in volumetric current density. In particular, note that the current distribution across the CL thickness is more uniform for the thin CL as seen previously. 72 Galvano-dynamic simulations.-For this set of simulations, a current profile shown in Fig. 3 is applied and the cell voltage is calculated. The profile is made up of the following steps: 0.2-1.0 A/cm 2 , 1.0-1.8 A/cm 2 , 1.8-1.0 A/cm 2 , and 1.0-0.2 A/cm 2 . The magnitude of the steps are chosen to be relatively high in order to excite the system dynamics. Similar to the voltage steps, these step changes are smooth and happen over a period of 1 second.
Overall, 16 simulations are conducted under the galvano-dynamic mode with variations in RH (60, and 90%), operating temperature (60, and 80°C), and cathode Pt loading (0.4, 0.2, 0.1, 0.05 mg/cm 2 ). The driest and coldest conditions used for the potentio-dynamic simulations could not be simulated in the galvano-dynamic mode with the selected current profile. This is due to the severe anode dry out with EOD that occurs during a step change in current density 12 and results in numerical issues under these dry conditions, where the membrane hydration is low prior to the increase in load.
The resulting voltage dynamics for all 16 simulations are shown in Fig. 10. The corresponding average water content in the membrane for all the cases are shown in Fig. 11 and the average liquid saturations in the cathode GDL are shown in Fig. 12. Similar to the analysis for the potentio-dynamic case, below we organize the discussion in terms of the current density step. Current density step from 0.2 to 1.0 A cm 2 .-After the first step in the current density, the voltage drops and a steady state value is achieved within 20 seconds. Only the case with T = 80°C and RH = 90% with a Pt loading of 0.4 mg/cm 2 demonstrates a slight undershoot that is a characteristic response of PEM fuel cells due to dry out of the membrane by EOD as mentioned earlier. 146 This dry out can also be observed in Fig. 11. It should be noted that mass transport limitations at higher loads can also contribute to this behavior. 152 However, membrane dry out is the only contributing factor in this case, since the operation is well within the ohmic region and mass transport effects are relatively insignificant. To finalize our discussion of the voltage response, it is worth mentioning that higher Pt loading is observed to consistently result in improved performance under all four simulated conditions at 1.0 A/cm 2 .
As for the GDL liquid saturation, only the coolest (T = 60°C) and most humidified (RH = 90%) condition results in vapor condensation after this step change in the current density (Fig. 12). This condensation continues throughout the entire 100 seconds hold at 1.0 A/cm 2 , where the higher Pt loadings result in slightly higher liquid accumulation in the GDL. This higher liquid saturation is due to the fact that the thicker catalyst results in higher cell voltage and reduced volumetric heat generation, which lowers the overall cell temperature by about 0.4°C. The average saturation of 0.2 agrees with the in-operando results by Banerjee et al. 153 They also propose fitting the time evolution of liquid saturation with a first order dynamic equation and obtain a time constant of 2.3 minutes for similar operating conditions using an SGL 25BC diffusion medium. Similar dynamics are observed by others as well. 154 Our results indicate a time constant of about 33 seconds, which is more than 4 times faster than that reported by Banerjee et al. 153 This discrepancy between the model predictions and experimental results may be attributed to differences in cell geometry, membrane and catalyst layers, and thermal properties assumed for SGL 24BC in this work. Current density step from 1.0 to 1.8 A/cm 2 .-The voltage response to the second step increase in current density is seen to be monotonically decreasing in most cases, while some cases (T = 60°C and RH = 60% in Fig. 10) exhibit the characteristic undershoot discussed earlier. The monotonic voltage decay is a sign of well humidified membrane and increasing mass transport limitations with slowly accumulating liquid water (e.g., T = 60°C and RH = 90% in Fig. 10 and Fig. 12). For the operating condition where a voltage undershoot is observed, the performance recovers to some extent after the step change with water back diffusion and rehydration of the anode side of the membrane. However, this voltage recovery is followed by a further decay as the cathode GDL floods with liquid water (Fig. 12).

RH = 90%
An important observation is the fact that at T = 60°C and RH = 90%, reducing the Pt loading from 0.4 to 0.2 mg/cm 2 seems to improve performance. This seemingly peculiar behavior is directly related to the changes in the membrane water content with CL thickness. In particular, Fig. 11 shows that at 60°C and 1.8 A/cm 2 , thinner catalysts consistently result in better membrane hydration. This is explained further when discussing the water balance in the cell later in this section. The higher membrane water content achieved at 0.2 mg/cm 2 reduces the ohmic drop, while the Pt reduction at this level does not impose significant mass transport issues. Therefore, the performance is improved. Yet another interesting observation is related to the impact of Pt loading on voltage dynamics at T = 80°C and RH = 90%. Specifically, we observe that immediately after the load increase, the highest Pt loading achieves the best performance. But this performance decays more rapidly than the cases with lower Pt loading, to the point that at the end of the 100 seconds hold at 1.8 A/cm 2 , the cell with a Pt loading of 0.2 mg/cm 2 has a higher voltage than that with a Pt loading of 0.4 mg/cm 2 . This is directly related to faster liquid accumulation in the GDL with a thick CL, which is due to lower cell temperatures (Fig. 12). The heat generation in the two thinnest CLs is high enough to inhibit any liquid accumulation (Fig.  12), which results in lower membrane water contents as seen in Fig.  11 and reduced performance.
Current density steps from 1.8 to 1.0 A/cm 2 and from 1.0 to 0.2 A/cm 2 .-For the simulated conditions, the step decreases in current density result in monotonic voltage increase with relatively fast dynamics. A hysteresis effect may be expected due to membrane hydration-dehydration and Pt oxide coverage dynamics. The former effect can be seen to some extent in the simulations results at T = 60°C and RH = 60%. The most notable feature of the voltage response to load decreases at this condition is a drop in performance after a considerable time delay. This is due to evaporation of liquid water reservoir in the cathode CL and the ensuing loss of membrane water (Fig. 11) and was discussed in detail for the potentio-dynamic simulations. Other conditions show almost no hysteresis, since they are well humidified. Furthermore, as mentioned in the Model Formulation section, the Pt oxide growth dynamics are neglected in our model, which further contributes to the lack of hysteresis in our results.
Further insight about the water balance in the cell can be gained by comparing the normalized membrane water fluxes (β) shown in Fig.  13. The figure illustrates β values (defined by Eq. 54) for different operating conditions and cathode Pt loadings at 119 seconds, which corresponds to the hold at 1.0 A/cm 2 just before the following step increase in the load. Some clear trends can be observed in these results. First, we note that increasing the RH increases β. This is due to the fact that the cathode CL typically has a high RH under all conditions due to the electrochemical water generation. A low RH dries out the anode CL and promotes water back diffusion, which explains lower β values. Another important observation is the impact of temperature on β. At lower RH conditions (RH = 60%), we observe that β increases with temperature. This is due to the fact that parts of the cathode CL remain subsaturated at low RH and high temperatures, which reduce water back diffusion and increases β. At higher RH conditions (RH = 90%), the cathode CL remains saturated at both low and high temperatures. Therefore, the water activity in the anode CL determines the driving force for water back flow. As this water activity diminishes at higher temperatures, the water back diffusion is more pronounced, which yields a lower value of β. Similar arguments can be used to explain the seemingly counterintuitive impact of Pt loading on the water balance. In particular, we note that at lower temperature, β increases with Pt loading and the resulting increase in CL thickness. This trend is reversed at higher temperature, where a higher Pt loading reduces β. The trend at T = 60°C can be explained by the fact that at this lower temperature, the CL is saturated for all Pt loadings. A lower Pt loading means that the water that is produced in the CL close to the MPL has a lower resistance to diffuse back to the anode, since the cathode CL is thinner and the diffusion path is shorter. This means that a thinner CL reduces β. When the temperature is increased, subsaturated conditions emerge in parts of the CL as mentioned before. Therefore, an increase in the temperature diminishes the driving force for water back diffusion and increases β. As thinner CLs generate more heat, the local CL temperature increases further, which in turn increases β. This result bears significance, as it shows that in addition to changes in local transport resistance, Pt loading impacts the performance by influencing the water balance in the cell. This overall observation partially confirms the hypothesis by Muzaffar et al., 77 who claimed that the performance changes with Pt loading reduction may mostly stem from a tipping water balance in the cell. However, their conclusion was based on the assumption that the CL is the main source of vaporization in the cell and a thinner CL makes the cell inherently more susceptible to flooding. On the other hand, the phase change rate is assumed to be relatively high in our model based on the experimental evidence in the literature that suggest the phase change kinetics should be fast enough not to impose any limitations. 143 This high rate of phase change allows the GDL to vaporize a relatively large amount of liquid. Nevertheless, our results also highlight the role of Pt loading and CL thickness in the cell water balance.
Our discussion of the galvano-dynamic simulations has focused on average response of the cell so far. However, the distribution of temperature, water, and reaction rates are also of critical importance. An example of such distributions is provided in Fig. 14 for high (0.4 mg/cm 2 ) and low (0.05 mg/cm 2 ) cathode Pt loading. The figure shows a rather significant in-plane temperature gradient in the membrane and F3173 CL, where a temperature difference of up to 4°C is observable. This temperature difference is due to the high current generation under the channel location along with the limited heat dissipation through the channel boundary.
Furthermore, we see that the cathode GDL has a very high saturation at this operating condition, which imposes significant mass transport limitations and results in limited current generation under the land area. It should be pointed out that a higher liquid saturation is observed under the land, while a rather uniform distribution is seen through the GDL thickness. Zenyuk et al. found a higher liquid saturation under the channel in a compressed GDL due to the in-plane porosity distribution resulting from land compression. 155 However, the temperature distribution in an operating fuel cell usually shifts the water accumulation toward the land region. 153 This flooding under the land has much less pronounced impact on the overall performance compared to flooding under the channel, 156 since the land region is already transport limited by the longer diffusion paths. In terms of the through-plane liquid distributions, there is evidence in the literature for increased saturation in higher porosity regions of the GDL. 113 Furthermore, Banerjee et al. 153 found the highest level of saturation to occur close to the land. Finally, capillary fingering is believed to be a major transport mechanism for liquid water. 155,157 The model in this work assumes a constant GDL porosity through the thickness, so the pooling effects cannot be captured. The rather uniform throughplane liquid distribution stems from the temperature distributions in our simulations. Finally, the macro-homogeneous model in this work does not allow for simulation of capillary fingering. Therefore, the liquid saturations predicted by the proposed model are only insightful on an aggregate level and detailed knowledge about the micro-structures are needed to obtain accurate distributions. 16 The non-uniform current generation pattern also affects water distribution in the ionomer phase. In particular, we see that under the land, the water content is more uniform across the thickness of the CCM, whereas significant gradients emerge under the channel, where current generation is high. Comparing the low and high Pt loading cases, we observe that the lower Pt loading results in slightly higher temperatures, which in turn reduce liquid saturation in the cathode GDL. Moreover, the current distribution through the thickness of the cathode CL is more uniform for the lower loaded CL as discussed earlier. 72 The preceding analysis provides some insight about the quasi steady state distributions of critical variables. To better understand the water transport transients during the following current step down, the distributions of liquid pressure in the cathode CL and the membrane water content before and after the step change are shown in Fig. 15 and Fig. 16, respectively, for the high Pt loaded CL (0.4 mg/cm 2 ). Fig.  15 shows that immediately before the step change at 219.5 seconds, liquid pressure is highest under the channel. The flow directions provided in the figure show an interesting pattern, where the liquid water is found to flow mostly toward the membrane in the land region and mostly toward the MPL in the channel region. This behavior is closely tied to a similar flow pattern in the membrane as seen in Fig. 16, where a recirculation is observed at 219.5 seconds. More specifically, it is seen that at the furthest location under the land, the membrane water flux is toward the anode. As we move closer to the channel location, the flux turns progressively toward the cathode.
As the load is reduced to 1.0 A/cm 2 , EOD is relieved and water back diffusion dominates during the transients (as was observed in Fig.  7 for the potentio-dynamic simulations). This results in the membrane water flux to be dominantly toward the anode (Fig. 16 at 221 seconds), which also pushes liquid water toward the membrane to compensate for the back diffusion ( Fig. 15 at 221 seconds). Two seconds after the step change, a smoother water profile is established across the membrane thickness and water back diffusion has diminished (Fig. 16 at 222.5 seconds). This creates a dominant flux toward the cathode. At the same time, a higher liquid pressure is observed in the cathode CL, where a stagnation front has emerged in the middle (Fig. 15 at 222.5 seconds). In particular, in the half of the CL close to the membrane, the flow is found to be toward the membrane, whereas in the other half to the right of the stagnation front, the flow is found to be toward the MPL. As time goes on and a quasi steady state is achieved at this reduced load, the stagnation front moves further toward the membrane (Fig.  15 at 230 seconds), while the membrane water flux turns toward the cathode throughout the membrane thickness ( Fig. 16 at 230 seconds). The existence of the stagnation front is in agreement with our earlier observation that a thick catalyst layer increases the resistance to water Figure 15. Liquid pressure (in [Pa]) and flow in the cathode CL during load decrease from 1.8 to 1 A/cm 2 at T = 60°C and RH = 90% (Pt loading of 0.4 mg/cm 2 ). From left to right: immediately before the load change, 0.5 seconds after the load change, 2 seconds after the load change, and the quasi steady state achieved after the load change. flow toward the anode. The transients during a step increase in the load are the reverse of those presented here.
The dynamics of in-plane current density distribution during load changes can also improve understanding of the transient phenomena.
To this end, the corresponding distributions for the hot and wet operating condition (T = 80°C and RH = 90%) during the load increase from 1.0 to 1.8 A/cm 2 are shown in Fig. 17 for both low (0.05 mg/cm 2 ) and high (0.4 mg/cm 2 ) cathode Pt loading. We see that immediately before the step change, both Pt loadings result in a relatively uniform in-plane current distribution. However, as the load is increased to 1.8 A/cm 2 the region under the channel tends to generate more current with both Pt loadings. It can be seen, however, that the in-plane distribution is more uniform for the lower Pt loading case. This is attributable to the fact that the thinner CL results in higher cell temperatures and lower Figure 17. In-Plane current density distribution for the galvano-dynamic simulations during load increase from 1 to 1.8 A/cm 2 at T = 80°C and RH = 90% for high and low cathode Pt loadings. From left to right: immediately before the load change, mid-way through the load change, immediately after the load change has completed, and the quasi steady state achieved after the load change.
) unless CC License in place (see abstract). ecsdl.org/site/terms_use address. Redistribution subject to ECS terms of use (see 207.241.231. 83 Downloaded on 2019-04-27 to IP RH = 60% RH = 90% Figure 18. In-Plane current density distribution difference between high-and low-loaded CLs during load increase from 1 to 1.8 A/cm 2 . liquid build up, which reduces the mass transport limitations under the land region. We also note that current generation under the land is low even immediately after the load increase and before any liquid accumulation. This agrees with experimental results by Schneider et al., who found the land region to be transport limited even prior to any liquid build up in the GDL. 158,159 To further investigate the impact of Pt loading on the dynamics of in-plane current distribution, the difference between the current densities with high and low Pt loadings during the same load change are shown in Fig. 18 for all simulated operating conditions. The figure shows that at 1.0 A/cm 2 the difference between the current distributions is relatively insignificant as was seen earlier. As the load increases, however, the cathode Pt loading and CL thickness seem to have varying impacts on the current distribution at different operat-ing conditions. In particular, under drier conditions (RH = 60%), the higher Pt loading and thicker CL result in improved performance under the land region immediately after the step change. At T = 60°C this performance enhancement under the land fades toward the new quasi equilibrium state, as more liquid builds up in the GDL when a thicker CL is used. However, at the higher temperature (T = 80°C) that inhibits liquid accumulation, the new equilibrium results in a distribution difference that is close to the one obtained immediately after the load change. When the RH is increased to 90%, the response is dominated by liquid accumulation dynamics in the GDL. More specifically, thicker CLs make the GDL more susceptible to rapid liquid build up, as they generate less heat. Therefore, the performance is diminished under the land region. This results in a particularly significant difference at T = 80°C and RH = 90%, where the thick CL prompts considerable GDL liquid saturation at 219 seconds, whereas the thin CL leads to a relatively dry GDL. Therefore, the high Pt loaded cathode CL shows a remarkable performance drop in the land region.

Summary and Conclusions
A comprehensive model that captures the most salient transient phenomena across the thickness of a unit cell is developed in this work. The model draws from and extends the existing models in the literature by incorporating state of the art reaction kinetics for the HOR and ORR, the mixed wettability model for porous layers, a consistent homogeneous model for the CL micro-structure, as well as the ionomer relaxation dynamics. The model predictions are compared with experimental data from the literature obtained through voltammetry and voltage step experiments under a variety of conditions and a good agreement is obtained.
The developed model is executed with different humidity and temperature conditions under both current and voltage control operational modes and varying Pt loadings in the cathode CL. The results of these simulations shed light into the transient processes that determine the dynamic response of PEM fuel cells to load changes. In particular, we have found the transient response to be dominated by water redistribution in the cell. The timescales of this redistribution are dependent on the operating conditions and are controlled by the membrane water uptake and two phase flow in the DM for dry and wet conditions, respectively. Furthermore, the modeling results suggest that changing the cathode Pt loading, and thereby, the cathode CL thickness, can influence the performance by affecting the water balance in the cell. Specifically, the thiner CL results in higher rates of heat generation on a volumetric basis while leading to a shorter diffusion path for water transport toward the anode. Our simulation results suggest that, based on the operating conditions, the combination of these effects lead to distinctly observable trends in normalized membrane water flux with respect to changes in the cathode Pt loading. Additionally, we have found that through its effect on water balance in the cell, the cathode Pt loading can have a profound impact on the transient response to load changes for some operating conditions. These findings can further improve understanding of the impacts of Pt reduction on various aspects of PEM fuel cell performance and its transient response. More broadly, the model can be used to develop further insight into spatiotemporal distribution of variables that are critical to performance and degradation.