Resolving Anodic Current and Temperature Distributions in a Polymer Electrolyte Membrane Water Electrolysis Cell Using a Pseudo-Two-Phase Computational Fluid Dynamics Model

Expanding upon our prior experimental work, we constructed a three-dimensional model of a polymer electrolyte membrane water electrolyzer using computational fluid dynamics. We applied the assumption of pseudo-two-phase flow, the flow of two phases with equal velocity. Experimental data were used to obtain parameters and to determine the conditions under which this model was valid. Anodic distributions of current density, temperature, liquid saturation, and relative humidity were obtained at various flow rates. The overall current density and temperature difference from inlet to outlet at the anode agreed strongly with experimental measurements under most circumstances. This verification allowed us to further examine the apparent gas coverage calculated from experimental and model temperature data. Results suggested a low liquid saturation and low relative humidity at the anode due to the consumption of liquid water and water vapor. However, we questioned the accuracy of the pseudo-two-phase assumption at low water feed rates. We concluded that the model was applicable to systems with liquid water feed rates greater than 0.6 ml min cm. Therefore, it is a fair screening method that can advise which operating conditions lead to excessive temperatures or drying at the anode, thereby promoting the longevity of the membrane and catalyst. © 2021 The Author(s). Published on behalf of The Electrochemical Society by IOP Publishing Limited. This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 License (CC BY, http://creativecommons.org/licenses/ by/4.0/), which permits unrestricted reuse of the work in any medium, provided the original work is properly cited. [DOI: 10.1149/ 1945-7111/abfe7b]


List of symbols
Computational fluid dynamics (CFD) is a versatile technique with numerous applications in the field of electrochemistry, especially flow batteries, fuel cells, and electrolyzers. Polymer electrolyte membrane water electrolysis (PEMWE) devices are becoming more widely commercialized as hydrogen consumers seek ways to reduce costs. The combination of electrolyzers and fuel cells provides an energy dense, flexible storage method as part of a diverse energy landscape. While the performance of PEMWE devices already make them competitive in a hydrogen market, 1 there is always a desire to optimize their design and operation.  Design optimization is most cost-effectively achieved by screening component geometries and material parameters with developed simulations, avoiding the fabrication of many physical prototypes. Computational methods can also assist the operation of existing devices. In research, modeling helps us to understand electrochemical systems. Before we can use computational methods for optimization, they must be validated by theory, experiments, and reproduction. The presence of the gas phase in electrolyzers reduces the effective conductivity of the electrolyte and the electroactive area of the catalyst. In alkaline electrolysis, this is because ions exist only in the liquid phase. CFD models, applied to electrolyte-fed verticalcolumn membraneless electrolyzers, considered two-phase flow and the impact of bubble curtains on effective electrolyte conductivity and electrode kinetics. Numerical two-phase flow profiles were computed by Hreiz et al. 2 and Le Bideau et al., 3 then scrutinized via comparison to experimental work. 4 El-Askary et al. 5 studied void fraction profiles, including the effect of bubble coverage on the electrode surface. The influence of bubbles on both effective conductivity and kinetics were included in an electrochemical model presented by Aldas et al., 6 which studied void fraction profiles. Twophase CFD models were developed also for alkaline diaphragm water electrolysis devices, which contain a porous separator that prevents gas mixing. Rodriguez and Amores 7 developed a validated model which included gas phase effects on effective electrolyte conductivity in the resistance overpotential. The electrochemical model by Mat et al. 8 additionally considered bubble coverage and ionic species transport. They reported the current distribution along the channel length. In PEMWE devices, the focus of this work, ions migrate through a solid phase, so while the effects of the gas phase are analogous, they are indirect. Ion transport and charge transfer are not confined to the liquid phase, but these mechanisms rely on sufficient humidification of the ionomer and catalyst surface. Therefore, it is important to quantify the liquid saturation and relative humidity at the anode in PEMWE devices. This motivated many efforts to model flow in PEMWE flow channels. Nie and Chen 9 used two-phase CFD with phase slip and drag forces to show how these phenomena could impact fluid velocities in electrolyzer channels, thereby impacting the oxygen volume fraction distribution. Olesen et al. 10 found little impact of phase interactions on flow distribution in their comparison of single-and two-phase models of flow through a high-pressure cell with interdigitated channels. Later on, Lafmejani et al. 11 demonstrated that a volume-of-fluid (VOF) model could be used to resolve different flow regimes within the interdigitated channels of the same high-pressure cell as Olesen et al. 10 While these authors studied the flow channels, Arbabi et al. 12 focused on modeling two-phase flow through porous media using VOF simulations of a lab-on-chip experiment. They were able to reproduce, with some accuracy, lab-on-chip experiments intended to represent porous transport layer (PTL) materials using a representative structure. These models demonstrated a capability to characterize twophase flow in PEMWE devices, but they did not resolve current distributions. Toghyani et al. 13,14 applied CFD to study high-temperature, gas-phase water electrolyzers, obtaining current distributions and fair agreement with experimental data. However, these models, as intended for their purposes, did not include water in the liquid phase, so effects of phase composition were not considered.
Two-phase flow in PEMWE cells was linked to the electrochemical reaction in a mathematical model developed by Han et al. 15 Capillary forces and viscous resistance governed flow while kinetic limitations were applied through a diffusion overpotential. This model accurately reproduced experimental observations of cell performance for different PTL thicknesses, porosities, and contact angles. However, being a one-dimensional, isothermal model, it was not intended to determine the current distribution or the local anode temperature.
Electrochemical models with two-phase flow in polymer electrolyte cells have been developed for water management in fuel cells. [16][17][18][19][20] Advanced three-dimensional (3D) fuel cell models include effects of capillary forces and porous media on two-phase transport. This is coupled through source terms or boundary conditions to a reaction, namely the oxygen reduction reaction, which occurs in the liquid and gas phases. However, similar work applied to PEMWE is needed. Herein, we propose a validated 3D CFD model that employs our present understanding of the two-phase anodic reaction in a liquid-fed PEMWE device.
Within the CL at the anode, the oxygen evolution reaction (OER) occurs at sites in proximity to both the catalyst and ionomer. At these reaction sites, liquid water is oxidized at the catalyst surface to form dissolved oxygen, protons, and electrons. Protons migrate to the ionomer and through the ionomer network toward the membrane, dragging water molecules with them. Electrons travel through the catalyst and support to the PTL. Near the catalyst surface, oxygen exists in a supersaturated state and so it nucleates to form a gas phase. Thereby, some of the catalyst surface is wetted, and some fraction is in contact with the gas phase. Upon the formation of a gas/liquid interface, there must be some evaporation if any liquid water is present. Water vapor diffuses through the oxygen gas toward the catalyst surface, where the OER occurs also in the gas phase. At sufficiently high current, the gas-phase reaction is limited by this diffusion, resulting in a concentration overpotential. This mechanism is challenging to deconvolute-the liquid saturation and relative humidity must be known-but this is a prerequisite for a model that can be used to engineer the PTL/CL interface where this mechanism takes place.
This work aims to couple pseudo-two-phase flow (PTPF) with a two-phase anodic electrochemical reaction to yield cell performance and temperature resembling that of an experimental cell with similar geometry. Pursuant to validated simulation results, we present the changes in steady-state anodic current density, temperature, liquid saturation, and relative humidity distributions with changes in water feed rate. Combining computations, our prior experiments, 21 and observations made by other authors, [22][23][24][25] this work also provides some speculation about the actual extent of liquid saturation and humidification at the anode. The data are used to suggest a minimum water feed rate conducive to reliable predictions of cell performance and temperature using the PTPF assumption.

Model Development
Geometry.-The geometry and polyhedral mesh are illustrated in Fig. 1. Figure 1a is a magnified image of the membrane and PTL mesh. Figure 1b adds the channel geometry. Figure 1c contains a cross-channel section of the polyhedral mesh inside the channels, PTLs, and membrane. The full geometry, including the plates and other structural components, is illustrated in Fig. 1d. The total number of computational cells in the entire geometry amounted to 3,643,729. A mesh independence study was performed on a simulation in wet anode/wet cathode mode with 80 ml min −1 water feed to each inlet at 2.0 V. The difference in current density between meshes with target surface sizes of half and double those of the original mesh was approximately 1.0%. The 2 cm × 2 cm cell consisted of 10 parallel channels, each with an approximate cross-sectional area of 1 mm 2 . There were 5 injection points perpendicular to the electrodes at the inlets to each flow field, and there were 5 identical exit channels at the outlets from each flow field. The PTL and membrane thicknesses were 1 mm and 0.2 mm, respectively. The latter dimension was based on the approximate thickness of a water-saturated Nafion® 117 membrane. 26,27 Computational efficiency was enhanced by using a small mesh size of about 0.0667 mm near electrodes and coarsening the mesh near the boundary separating the channels and plates. The simplified model geometry did not include rubber gaskets, stainless steel fittings, or water feed and exhaust lines, the aim being to use an efficient computational mesh that leads to fair performance and temperature predictions. and energy are provided in Eqs. 1-20. Refer to the nomenclature section for the definitions of symbols, subscripts, and superscripts. Parameters are provided in Table I and additional relationships are provided in Table II. The eleven regions are labeled by Roman numeral in Fig. 1. These numerals will be used when appropriate to specify regions and interfaces.
In the following equations, superficial variables are signified with a subscript S and are defined as φ S = χφ, with χ being porosity and φ being any property. Volume-weighted average properties such as the density ρ and viscosity μ are defined as , l g , å j e j = P= P P in which ∏ is either liquid or gas and ε is the volume fraction. The built-in steady-state equations for the conservation of mass, momentum, and species were applied to the channels and PTLs. e being the relative permeabilities of the liquid and gas phases, respectively. The capillary pressure is expressed below using a Leverett function J applicable to hydrophilic porous media (see Table II), which can be found elsewhere. 19 ¹ to allow phase change between liquid water and water vapor. Note the use of Fick's Law in the first term on the righthand side, which is an approximation. 42 A no-slip condition was assumed at the fluid-solid interfaces. At the channel inlets, the velocity of liquid water feed was specified based on the volumetric flow rate Q and the cross-sectional area of the inlet A {in} : At the electrodes, species were generated and consumed. Mass fluxes due to the electrochemical reaction, electroosmotic drag, and phase change are given in Eqs. 9, and 11, respectively, with ν, M, F, i , ̲ k, ̲ nd, , ̲ x and y respectively representing the stoichiometric coefficient, molecular weight, Faraday's constant, current density, the through-plane unit vector, the electroosmotic drag coefficient, the inner normal unit vector, and mole fraction. Equation 11 assumes that evaporation in the CL was sufficient to maintain the mole ratio of water vapor. The simplifying assumption for PTPF is that both phases have equal velocity. This led to one governing equation each for mass and momentum conservation, easing the computational load, but an effective phase diffusivity must be defined, Anodic equilibrium potential of the gas-phase reaction Anodic equilibrium potential of the liquid-phase reaction e e e -+ 19, 20, 33

Weight fraction of water vapor at saturation
The capillary diffusivity was defined in Case 2 of Eq. 12 for Regions III and IV, slightly modified to have a lower bound for numerical stability. 43 Case 1 of Eq. 12 is the phase diffusivity in the channels. The volume-weighted average Schmidt number was determined as follows: Ohm's Law and conservation of charge were applied to the membrane, plates, and PTL regions, with σ denoting conductivity.
The boundary conditions for the potential E at the plate contacts and anode were as follows: The area-specific thermodynamic and kinetic resistance at the anode interface, R tk,{an} , was adjusted to control the current density using the Nernst equation and a modified Tafel equation, which is explained in the next section. The conservation of energy was applied to all regions, some having region-specific terms, in which ĥ is mass-specific enthalpy of formation, k is thermal conductivity, and T is temperature. Equation 17 includes terms for convection, viscous stress heating, and an enthalpy source term S I IV for evaporative cooling in regions I-IV. Ohmic heating is considered in the membrane while it is assumed negligible elsewhere. All other regions have only the first term on the right-hand side, associated with conduction, equal to 0. The following heat flux conditions were enforced at the inlets, outlets, and electrode boundaries: I II  S  I II  I II  in  , , , ,

I II in
The right-hand side of Eq. 19 includes, in order, a boundary source term for reaction heating, evaporative cooling, heat transfer through the membrane with water, and conduction through the membrane. At all outer surfaces besides the channel inlets and outlets, a convective cooling boundary condition was applied with U as the convective heat transfer coefficient: In the wet anode/dry cathode configuration, zero-flux conditions were applied to the cathode inlet for all relevant physics.
Electrochemical equations.-As shown in Fig. 1, the electrodes were assumed to have infinitesimal thickness and were defined at interfaces. At the cathode interface, the hydrogen evolution reaction (HER) occurred: It was assumed that the HER was ideally non-polarizable and that the equilibrium potential was fixed at 0.0 V, resulting in no potential drop at the cathode interface. Meanwhile, at the anode interface, the OER occurred in two phases, each with different standard equilibrium reduction potentials.
The Nernst equation was applied at the anode to determine the equilibrium reduction potential, in which  is the gas constant and a i,Π is the activity. The activity of water vapor was equal to the relative humidity h r , the activity of oxygen was equal to its gas-phase partial pressure, and the liquid water activity was set to 1. Once the equilibrium potential was determined, it was possible to solve for the kinetic and resistance overpotentials. R tk,{an} , in Eq. 16, was determined at each face on the boundary iteratively using a modified Newton-Raphson method. The total through-plane current density at the anode interface i {an} was calculated using two Tafel approximations, one for each phase, with that of the gas phase being concentration-dependent: 44 In the above equations, i 0 is the exchange current density 21 and α is the transfer coefficient, the latter of which has been found to apparently vary linearly with temperature. 21  To conserve mass, the phase change rate cannot exceed the flow rate of liquid water or the maximum uptake rate of the gas phase, that is: assuming that the weight fraction of the liquid phase w l @ 1 and that ρ g is uniform in the control volume. The phase change factor β was adjusted to match performance between the simulation and experiment, which allowed us to suggest relative humidity values at the anode. While β was used as a fitting parameter, it does have physical significance, as it should increase with increasing specific interfacial area between the liquid and gas phases. Model integrity was ensured by using the same β value in all simulation runs. A β value of 0.07 was used in porous media, which have high phase interaction areas. β was set to 0.001 within the channels. The sensitivity of β within the PTLs, especially at low water feed rates, highlights the importance of quantifying phase interaction areas in electrochemical porous media.
Numerical procedure.-Optimal solutions for this system were obtained using Siemens Simcenter Star-CCM + 15.02.009-R8, which utilizes the finite volume method. An overview of the general computational sequence is provided in Fig. 2. The initial condition for potential was set using Eqs. 21-28 to estimate the potential on the membrane side of the anode interface. A linear potential profile was assumed within the membrane initially. The convergence scheme began with obtaining a solution for the potential, then activating the other solvers. First, current was calculated using the initial condition for potential, then the potential solver was allowed to converge on its own. The final under-relaxation factors (URFs) for the potential, velocity, pressure, species, fluid energy, and solid energy solvers were set to 1.0, 0.7, 0.3, 0.995, 0.995, and 0.9995, respectively. Upon activating the flow, species, and energy solvers, flow solver URFs were set to their final values while species and energy URFs were ramped from 0. After the ramp, the gas diffusion limitation due to the presence of liquid, then flow compressibility, were activated sequentially. Convergence was said to be achieved once the rate of change in current density was less than 10 A m −2 per 1000 iterations and slowing, which occurred long after the residuals stopped decreasing. At low flow rates, species and fluid energy URFs were reduced in response to numerical instability caused by anode evaporation-see Eq. 11-overcompensating for large fluctuations in relative humidity, leading to run-away oscillation.
Model interpretation.-Potential components were determined using the same assumptions as those for the experiment; only the liquid-phase equilibrium potential was compared to the experimental equilibrium potential. The resistance overpotential was determined via the following equation: with A being the cross-sectional area of 4 × 10 −4 m 2 and V being the volume of the membrane. The nominal activation overpotential was then calculated as the difference between the cell potential and the IR-corrected potential. During experiments, 21 the CL temperature was unmeasured. It was cautiously assumed that the reaction temperature was that of the anode outlet. We used this temperature to calculate an apparent gas coverage that helped us to interpret the experimental results. It was acknowledged that the actual anode temperature may have been higher than the measured outlet temperature, which was a possible explanation for the observed decrease in apparent gas coverage with a reduction in stoichiometric flow rate. We defined the apparent gas coverage from simulations in the same manner as in Eqs. 12 and 13 of our prior experimental analysis, 21 is the anodic overpotential in the liquid phase adjacent to a water-saturated gas phase. Succinctly, Θ app is the fraction of the exchange current density that is apparently deducted due to the presence of gas, assuming that uncompensated CL ionomer resistance is negligible. The actual gas coverage from simulations was defined in this work as the volume fraction of the continuum that cannot be consumed: In the results and discussion, Θ app values from our previous experiments are corrected based on anode temperatures computed in the simulations. Θ app and Θ act are then used to explain the simulation results at low flow rates.

Results and Discussion
In this section, we validate the 3D model using experiments and then analyze the results. The ideal model requires knowledge of three variables at the anode: temperature, ε l , and h r . High temperatures enhance kinetics and reduce membrane resistance while slightly reducing thermodynamic favorability. ε l is related to the fraction of catalyst that will facilitate the splitting of liquid water. h r determines whether the gas phase enhances performance due to the relatively low equilibrium potential of the gas-phase reaction (high h r ) or reduces performance due to mass transfer limitations (low h r ). In the study of these three variables, we first show the range of validity of the isothermal assumption at flow rates of 80 ml min −1 , which supports the empirical relationships in Eqs. 23 and 24. The distributions of current density, temperature, liquid saturation (ε l ), and relative humidity (h r ) at various flow rates are then provided and detailed. Next, we compare the measured temperatures to the model temperatures at various flow rates. An accurate thermal model allows the use of anode temperatures suggested by the simulations and Eq. 34 to correct the Θ app calculated from experimental data. Then, the suggested values of the ε l and h r from simulations are assessed and we determine the accuracy of the PTPF flow assumption. Finally, we use this evaluation to reach conclusions and identify subjects of future study. Refer to Table III for the properties of the  selected anode PTLs, PTL1 and PTL2, and use Table IV to reference Series A (wet anode/wet cathode, 80 ml min −1 per half cell, PTL1, variable potential) and Series B (wet anode/dry cathode, 1.9 V, variable feed rate, PTL1 & PTL2) of the operating conditions used for validation and in the figures.
Isothermal and non-isothermal cell performance.- Figure 3 compares the cell performance and potential components measured in experiments with the simulated performance, both with and without the isothermal assumption, under operating condition Series A. The non-isothermal simulation (NITS) gave a more accurate prediction of current density at 2.0 V, at which temperature effects were greatest. Decomposing the potential, we found that its components were well predicted by the models as well. The NITS became more important for performance estimations when the cell was operated under conditions that induced temperature change-Series B-as presented later on. The current densities computed by the isothermal simulation (ITS) deviated from those of the NITS by less than 1.5% at potentials up to 1.8 V. At 1.9 and 2.0 V, the deviation increased to 2.4% and 3.7%, respectively, with the ITS slightly underpredicting current density. The small difference between the two models under these operating conditions suggests that the extraction of kinetics from experimental data, which assumed temperature to be independent of potential and current,  was valid. Therefore, Eqs. 23 and 24 can be used to predict local temperature effects. As fluid flowed through the channels from inlet to outlet, it was heated by the anodic reaction and ohmic heating in the membrane. This is shown in Fig. 4, which shows the temperature of the fluid in the anode channels, including the inlet and outlet manifolds. Under Series A5 conditions (wet/wet, 2.0 V, 80 mlml min −1 per half cell, PTL1), despite the relatively high potential, the fluid temperature rose minimally, with the highest temperature near the end of the channel. This temperature change led to nonuniformities in other variables.
Distributions of current density, ε l , and h r from the ITS and NITS in Series A4 (wet/wet, 1.9 V, 80 ml min −1 per half cell, PTL1) are juxtaposed in Fig. 5. In the current distributions, the flow field pattern is clearly visible, with the current density being highest near the land area of the flow field as opposed to the channel, where current density is lowest. This is due to geometry-dependent pressure variations, with the highest pressure existing where fluid flows toward the anode from the channel. High pressure leads to lower electrolysis efficiency, and thus lower current density at constant potential. ε l decreased along the channel from inlet to outlet, which is due to gas generation and liquid water consumption. h r also decreased from inlet to outlet due to oxygen generation and water vapor consumption.
The current density increased as expected in the NITS due to temperature rise. There is also a localized temperature effect in the NITS leading to increased current density near the outlets. This is caused by the increase in temperature along the channel, which was demonstrated in Fig. 4. ε l was slightly lower in the NITS due to gas phase expansion and more gas generation at higher temperatures. h r changed insignificantly between models.
The pressure variations at the anode can be explained through flow visualization across the channels. Figure 6a shows the gas volume fraction (ε g ) distributions obtained from Series A1 (1.6 V) and A4 (1.9 V) across the channels with vectors indicating the direction of fluid flow tangent to the cross section. These vectors point toward the anode near the middle of the channel and point away from the anode near the ribs of the flow field and the channel edges. There were opposing forces acting on fluid in the PTL: pressure and capillary forces. Fluid flowed toward the anode where   the capillary force from Eq. 4 exceeded the total pressure gradient and flowed away from the anode otherwise. The outcome of this was that fluid generally flowed from the middle of the channels to the anode and from the anode to the edges of the channels. This was expected behavior based on experimental observations by Leonard et al., 47 who visualized the locations of gas bubbles in the PTL with respect to the flow field geometry and quantified the velocities of the liquid and gas phases. The rather small gradient in volume fraction within the PTL was unexpected because it had been previously observed that the liquid saturation decreases substantially in the through plane direction from the channel to the anode and is close to 1 near the flow field. [22][23][24][25] The magnitude of the phase diffusivity from Eq. 12, which reached values ∼10 −4 m 2 s −1 throughout the PTL, is the likely cause. As shown in Fig. 6b, ε g increased along the length of the channel and increased with cell potential as expected due to gas generation.
Comparison of anode feed rate sensitivity between experiments and simulations.-We used the NITS to predict cell performance under Series B conditions. Figure 7 shows experimental and simulated current densities at various flow rates. The model reproduced the increase in current density with reduced flow rate that was observed in our prior experiments and elsewhere. 48 At water feed rates of 2.4 ml min −1 and higher, the computed values agreed with experiments within 4% error. However, the experimentally observed optimal water feed rate was not predicted by the model, indicating that the predicted ε l or h r may be too low. The model was also somewhat sensitive to bulk properties of the anode PTL. Changing from PTL1 to PTL2 led to a 300-700 mA cm −2 decrease in current density that was not observed experimentally. These discrepancies may be attributed to the PTPF continuum model used, the assumption that PTL properties were uniform and isotropic, or the evaporation model.
In Fig. 8, current density, temperature, ε l , and h r distributions at the anode in Series B1-a (80 ml min −1 , PTL1), B3-a (5.0 ml min −1 , PTL1), and B5-a (0.82 ml min −1 , PTL1) are compared. The current density increased upon decreasing the feed rate from 80 ml min −1 to 5 ml min −1 due to temperature rise, but there was a shift in the location of highest current density, which was attributed to the reduced convective heat transfer. At low flow rates, the highest temperature was centralized because the convective component of heat transfer, which carries heat along the channel, was not as prominent. From 5 to 0.82 ml min −1 , the current density decreased overall and became less uniform. It was highest near the inlets and lowest near the outlets. This is explained by the reduced ε l and h r , expected symptoms of insufficient liquid water feed. The ε l and h r decreased along the channel as expected at all flow rates. These variables decreased and became less uniform at lower water feed rates as more water was consumed relative to the amount supplied.
Inferences made from gas coverage at the anode.-In experiments, temperature was measured in the space surrounding the cell and at the inlets and outlets. The average temperatures at the inlet and outlet boundaries in the simulations were recorded and the difference was compared to experimental results from Series B in Fig. 9. The experimental temperature change was defined as the difference between the inlet and outlet temperatures at 80 ml min −1 and the difference between the oven and outlet temperatures at lower feed rates. The rationale was that the inlet temperature in the simulations was fixed at 328.15 K, but we observed an increase in this temperature at lower feed rates in experiments due to reduced convection, allowing conduction from the cell through the inlet line to influence the inlet temperature. At feed rates of 10 ml min −1 and below, the inlet temperature in experiments was controlled by the oven, which was held at the cell temperature setpoint. Upon comparison of the temperature change, there was excellent agreement with error of less than 0.4 K at feed rates of 2.4 ml min −1 and higher. At 0.82 ml min −1 , the experimental temperature change was much less than predicted, but this was likely caused by conductive heat transfer through the outlet line in the experiment. The agreement in temperature change supports the anode temperatures given in Fig. 8 and merits the use of the simulation anode temperatures in Eq. 34 to calculate the Figure 5. Current density, gas volume fraction, and relative humidity distributions at the anode interface (visible in the top extension as a yellow highlight) for isothermal and non-isothermal models in Series A4 (wet/wet, 80 ml min −1 per half cell, PTL1, 1.9 V). Approximate locations of inlets and outlets, respectively, are indicated by ⊗ and ○. experimental Θ app . Figure 9 also compares the experimental temperatures used in Eq. 34 to the average anode temperatures from simulations, the latter of which were much higher. Previously, we used a preliminary assumption that the anode plate temperature near the outlet was nearly equal to the anode temperature. Simulation results showed that this was certainly false. Figure 10 compares the gas coverage calculated from experimental data from Series B, the experimental gas coverage corrected using anode temperatures from the simulations, and the gas coverage calculated from Eq. 16. It is likely that the observed decrease in Θ app with reduced flow rate was caused by the anode temperature assumption. When the experimental gas coverage was corrected using anode temperatures from the simulation, this drop in gas coverage was eliminated. The Θ app from the simulations was much less than the corrected experimental values at moderate feed rates, which was likely due to a gas phase equilibrium potential that was lower than that of the liquid phase. There is a gap in Θ app between the two anode PTLs in the experiment and simulation. This appears to be consistent with experiments but note that the gap seen experimentally was attributed to the larger surface grain size of PTL2 hindering diffusion while the gap in simulations was due to the lower permeability of PTL2. Based on experimental findings, 21,49 bulk permeability was not considered a significant factor in performance at the magnitudes of interest. In this CFD model, the surface grain structure of the PTL was not considered, so we did not expect to see much variation between PTL1 and PTL2.   The increase in Θ app with decreasing stoichiometric feed rate starting at about 100 is due to an overpredicted optimal feed rate. ε l and/or h r may be too low at the anode interface at low flow rates. At 1.9 V, the average ε l ranged from as high as 0.38 at 80 ml min −1 to as low as 0.04 at 0.82 ml min −1 . The high value is lower than values measured by other authors. Satjaritanun et al. 22 used operando X-ray computed tomography (CT) experiments to measure ε g values between 0.2 and 0.4 using a 0.320 mm activated carbon PTL. Thinner anode PTL materials such as the 0.19 mm Toray paper used by Lee et al. 23 led to a consistent steady-state ε g of 0.2 at the anode. Seweryn et al. 24 observed ε g between 0.4 and 0.5 at the anode via neutron radiography (NR) when using a 1.2 mm titanium PTL. Zlobinski et al. 25 observed even higher steady-state ε of up to 0.6 in a 1.0 mm PTL at 2 A cm −2 with NR.
Overpotential isotherms representing h r {ε l } when Θ app = Θ act , derived from experimental data from Series B and Eqs. 34 and 35, are plotted in Fig. 11. h r decreases with increasing ε l at constant overpotential because if one variable is low, the other must be high to have similar overpotential. No trend was observed in the isotherms with stoichiometric feed rate, which is given in the Fig. 11 legend. Θ app can be greater or less than Θ act ; if greater, it indicates that the CL ionomer resistance is non-negligible. If less, the overall equilibrium potential, Θ act and its relationship with Θ app is unknown experimentally, but if they were equal and experimental ε l values were within the literature range, then h r is likely in the range from 0.2 to 0.7, indicated by where the isotherms pass through the emerald-colored range of literature ε l values from NR. Average h r values computed in simulations ranged from 0.03 to 0.13. If the computed h r is accurate, then either the computed ε l is much lower than it should be or there is a significant effect of E an g eq , { } on performance. This assessment assumes that the CL ionomer resistance is negligible and ε l in the CL should be equal to the NR-observed ε l . Considering the very low feed rates and the heterogeneity of the PTL/CL interface, this might not be the case, especially in parts of the CL sandwiched between the PTL grains and the membrane. This underscores the need for accurate modeling of two-phase flow throughout both the anode PTL and the CL.
Based on simulation results and experimental validation, we found that the PTPF model is valid for predicting current and temperature distributions at feed rates higher than 2.4 ml min −1 for a 4 cm −2 cell, or 0.6 ml min −1 cm −2 , which is practical considering the use of higher feed rates to control stack temperatures. 50 Future developments.-The objectives of future work will be to improve the through-plane liquid saturation profiles in the PTL to match experiments, which may lead to lower optimal feed rates and thus more accurate current distributions. A phase change model that is built on fundamentals rather than parameters may be considered. Humidity effects on ionomer and membrane resistance are planned for future work. We intend to develop a model incorporating effects of the PTL structure on two-phase flow, resistance, and catalyst utilization.

Conclusions
An experimentally validated 3D CFD model was used to determine current and temperature distributions in a PEMWE device. The PTPF model documented in this work was a computationally efficient method of determining cell performance at high feed rates, making it a feasible design tool for temperature management in large stacks. Overall cell performance and temperature were determined with low error at flow rates greater than 0.6 ml min −1 cm −2 . The model's ability to resolve current and temperature distributions allows us to pre-determine the temperature of the anode before operating a cell. It thereby possesses utility in temperature management applications, helping to identify operating conditions that lead to unfavorably high local temperatures. The anodic ε l and h r values suggested by the simulations were plausible, but possibly underestimated, so the model may also be used to conservatively predict drying conditions. This work motivates the pursuit of models that can determine ε l and h r with more confidence and include the uncompensated CL ionomer resistance. Figure 11. Anode overpotential isotherms for low flow rate experiments using PTL1 (blue) and PTL2 (red), obtained by assuming Θ app = Θ act . Values in the legend are the experimental stoichiometric liquid water feed rates. Yellow shading highlights simulation ranges for relative humidity and liquid saturation while green shading highlights the range of liquid saturation observed in literature. If actual ε l and h r values exist in the region below the isotherms, it is due to the low equilibrium potential of the gas phase in simulations. Above the isotherms, high CL ionomer resistance is indicated.