Computationally efficient pseudo-2D non-isothermal modeling of polymer electrolyte membrane fuel cells with two-phase phenomena

In this article, a computationally efficient pseudo-2D model for real-time dynamic simulations of polymer electrolyte membrane fuel cells (PEMFCs) is developed with a specific focus on water and thermal management. The model accounts for temperature dynamics, two-phase flow and flooding in the diffusion media, and membrane water crossover as well as absorption and desorption processes. Computational efficiency is achieved by leveraging the disparate time scales within the system dynamics, in addition to exploiting the large aspect ratio of the cell layers to create a spatio-temporal decoupling. Taking advantage of such decoupling, the model yields a computationally efficient solution while providing detailed information about the state of water and temperature throughout the cell. Through this approach, the current implementation of the model is found to be about twice faster than real time. Moreover, a case study is carried out where different mechanisms contributing to overall water balance in the cell are investigated. The results are shown to be in qualitative agreement with published experimental data, thereby providing a preliminary validation of the modeling approach. Finally, using the modeling results, an equivalent electrical circuit model is proposed to help elucidate water transport inside various cell layers. © The Author(s) 2016. Published by ECS. 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/2.0871613jes] All rights reserved.

Real-time estimation, prediction, and control of cell hydration and temperature distribution is essential for optimizing the performance of polymer electrolyte membrane fuel cells (PEMFCs), as well as avoiding critical conditions and mitigating cell degradation. These applications necessitate mathematical models that not only run in real time, but also incorporate the important physical phenomena related to water transport and thermal management. However, including such phenomena comes at the cost of higher computational requirements, resulting in a trade-off between model accuracy and computational speed, which must be carefully balanced based on the desired application. As a result of these competing requirements, developing mathematical models that achieve a balance between the needs for high fidelity and low computational demand remains an active area of research.
Within the context of this paper, a fuel cell model is considered to have high fidelity if it incorporates the following phenomena: i) 3-D effects including anisotropic material properties 1-3 to resolve transport phenomena in all physical directions; ii) transient behavior; iii) detailed multistep hydrogen oxidation reaction (HOR) and oxygen reduction reaction (ORR) kinetics; 4,5 iv) multiphase flow in gas channels and porous media; v) non-isothermal effects; 6 and vi) multicomponent diffusion. 7,8 A more detailed explanation of these considerations follows.
In terms of dimensionality, 3-D models are of highest fidelity, because they are capable of capturing transport in both through-themembrane and along-the-channel directions and also account for the channel-land effects in the third dimension. Moreover, these models can easily incorporate the anisotropic material properties of various cell layers, which is especially important in the case of gas diffusion layers (GDLs). 9 Numerous such models have been reported in the literature, which usually utilize commercial computational fluid dynamics (CFD) or finite element method (FEM) packages. [10][11][12][13] On the other hand, these models are computationally very expensive and do not lend themselves to real-time applications. To address this issue, many models have been proposed, where the aspect ratio of cell layers are taken into account to decouple one dimension from the other(s). The so-called pseudo 2-D 14 and pseudo 3-D models 15,16 are obtained with this approach. Kim et al. 17 conducted numerical experiments, where they showed a pseudo 3-D model achieves the best balance between computational efficiency and model fidelity. Nevertheless, significant computational reductions can be achieved with pseudo 2-D models, where the channel concentrations are found and used to serve as spatially varying boundary conditions for the through-themembrane sandwich model. This approach is specifically useful for straight-channel flow fields, where the usual corner effects in interdigitated and serpentine flow fields are not present. Moreover, pseudo 2-D models may be calibrated to account for the channel-land effects through the use of effective resistances at the channel interface with the diffusion media. 17 In terms of temporal behavior of the fuel cell system, transient models are required to capture various dynamic phenomena in the cell that occur on multiple time scales and have a profound effect on its performance. Order of magnitude analysis has been reported in the literature to show the time constants associated with various dynamics in the cell. Mueller et al. 15 presented such an analysis, where they found the time constants for gas diffusion in porous media and reaction kinetics to be on the order of a few milliseconds. They also associated a time constant of 2 × 10 −7 s to electrochemical double layer charge and discharge. A more detailed analysis based on model order reduction (MOR) techniques was presented by McCain et al., 18,19 where they also found the dynamics associated with gas species diffusion to be very fast compared to other transients. These studies suggest that fast dynamics of gas transport in the fuel cell can be neglected with relatively small errors when time scales of larger than 10 milliseconds are of interest.
In terms of reaction kinetics, the Butler-Volmer (BV) equation is traditionally used to describe activation overpotentials on both anode and cathode electrodes. 20 However, recent studies have shown that the actual reaction kinetics in fuel cells do not follow BV model accurately 21,22 and therefore multistep models are better fitted to experimental data. 4,5 Water and thermal management of the fuel cell are very important aspects that are incorporated into many of the models in the literature. Temperature and two-phase flow effects are highly coupled through temperature dependence of saturation pressure and gas transport as shown by Weber et al. 23 This was experimentally validated in neutron imaging studies of Hickner et al., [24][25][26] where they found less liquid water residing in the cell at higher temperatures. In general, Journal of The Electrochemical Society, 163 (13) F1412-F1432 (2016) F1413 increasing the temperature results in a lower open circuit voltage and lower relative humidity in the cell as was confirmed in the modeling study of Bhaiya et al. 6 Furthermore, temperature gradient across various cell layers can significantly alter water distribution inside the fuel cell (through effects such as thermos-osmosis in the membrane 1 and phase-change-induced flow 27,28 ), thereby affecting its performance. Moreover, multicomponent diffusion can become an important factor at some operating conditions. Hickner et al. 25 found a loss of performance in the region where less liquid water residue was found in the cell at high current densities. This observation was attributed to additional resistance to oxygen transport due to high flow rates of water vapor from the catalyst layer towards the channel. These results indicate the importance of multicomponent transport for some operating conditions. However, for most of the operating range, binary mixture approximations are valid and Fick's law of diffusion can be applied with small errors.
As mentioned earlier, including all these effects results in increased computational costs, whereas lower costs are desired for real-time applications. To this end, many computationally efficient models have been proposed in the literature. Pukrushpan et al. 29 developed a 0-D isothermal and dynamic model for system-level control tasks. The model is highly efficient, but it does not provide any spatial resolution of the variables. The simplistic nature of the model limits its predictive capability to a narrow range of operating conditions. Furthermore, it requires extensive parameter identification to fit the modeling results to experimental data. Nevertheless, it remains one of the main models that is used in the control community. Philips et al. 30 proposed a reduced model for nonlinear model predictive control applications, where they used representative elementary volumes (REVs) to reduce partial differential equations (PDEs) into ordinary differential equations (ODEs). The model was claimed to simulate a 30 min time window in less than a second. However, lack of two-phase flow consideration caused the model predictions to deviate from experimental results at medium to high current density values. Golbert et al. 31 presented a pseudo 2-D model of PEMFC with emphasis on reactant flow rates in the channels. The only dynamical state of the model was the cell temperature. The model did not account for two-phase flow and other transient behavior and transport in the through-the-membrane direction was not resolved spatially, resulting in unrealistic current density distributions. The initial model was further reduced using continuously stirred tank approach to yield a simplified version of the model suitable for adaptive control applications. This reduction resulted in further decrease in model accuracy, thereby affecting its predictive capability.
There have also been numerous efforts at bridging the gap between high fidelity and computationally efficient models. Grotsch et al. 32 developed a reduced order 1-D model from their earlier modeling study. 33 The model accounted for two-phase effects in the diffusion media and catalyst layers and the dual pore size distribution in the membrane that was first proposed by Weber. 34 However, the isothermal assumption prevented the model from capturing temperature distribution and its effects on water transport in the cell. Moreover, the model did not account for reactant consumption along the flow channels. In other similar efforts, McCain 19 and Siegel 35 developed 1-D models that account for two-phase flow in the diffusion media as well as membrane hydration. Similar to previous studies, these models were isothermal and did not account for reactant consumption along the channels. Another study by Promislow et al. 36 presented a pseudo 2-D dynamic model that included the reactant channels and the cell sandwich, where the only dynamic states were those related to liquid water in the cell as well as the two-phase front location. Based on the profiles of liquid water in the diffusion media, they proposed tracking of the two-phase front instead of actual liquid accumulation in the diffusion media. The main drawback of this approach was its inability to capture water redistribution due to temperature gradients in the cell.
Thus, to the best of the authors' knowledge, no model is reported in the literature that simultaneously accounts for transient behaviors of the fuel cell system, two-phase effects and their coupling with spatial temperature distribution, and reactant consumption along the flow channels in a computationally efficient manner. The current work aims at filling this gap. Its goal is to include the abovementioned phenomena in a faster-than-real-time model that can be used for online estimation to inhibit cell degradation and prolong the life of the fuel cell system. This is done by taking advantage of spatial decoupling, disparate time scales and efficient numerical algorithms.
The rest of the paper is organized as follows: First, the mathematical model in different cell domains is developed in detail. Numerical implementation section presents the required boundary conditions as well as the decoupling approach adopted herein to achieve computational efficiency. A discussion of the numerical schemes used in the model is also provided. Next, some modeling results are presented followed by the discussion of the results and comparison to existing literature. Summary and conclusions are given in the final section.

Mathematical Model Development
The fuel cell model presented here is a semi-empirical, dynamic, pseudo-2D cell-level model that is developed for a straight channel configuration with co-flow reactant streams. The model consists of a submodel for flow channels, a cell sandwich model that describes reactant and water transport across various cell layers, and temperature and voltage submodels. Within the pseudo-2D domain, weak spatial coupling is utilized between the along-the-channel and through-themembrane directions. The only coupling between these two directions is achieved through boundary conditions. Fig. 1 shows the modeling domain for this study, which consists of the polymer electrolyte membrane (PEM), anode and cathode GDLs, and gas channels (GCs). The microporous layer (MPL) is not included in this model. MPLs are shown to affect the water balance in the cell and reduce contact resistances. 37  by which the MPL improves performance. 40 Moreover, interfacial phenomena (e.g. cracks in the MPL surface) at the MPL/GDL boundary have resulted in discrepancies between experimental and modeling results for detailed MPL effects on water transport. 41 Due to these difficulties, the MPL is not explicitly considered in this model, and it is assumed that its effects can be captured through effective material properties of the GDL. Additionally, as seen in Fig. 1, the regions under the lands are not considered in the modeling domain. The land effects can be captured through the use of effective resistances at the channel/GDL interface. 17 Furthermore, the 2-dimensional transport inside the anisotropic diffusion media can be modeled by a 1dimensional transport in an isotropic media using conformal mapping to find an effective diffusion thickness. [42][43][44] Finally, the land effects may be captured by solving the through-the-membrane direction twice with different boundary conditions as was suggested by Zaglio et al. 45 For the purpose of this study, it is assumed that the land effects can be captured with the effective convective resistance at the channel/GDL interface.
Other assumptions used for the model are as follows: r Catalyst layers are assumed to act as thin interfaces where electrochemical reactions occur (i.e. reaction distribution is assumed to be uniform).
r All gases are assumed to behave ideally. The range of operating conditions in a PEM fuel cell system permits this assumption.
r Total gas pressure in the through-the-membrane direction is assumed to remain constant.
r Due to very high interstitial area between solid and fluid phases, thermal equilibrium between the two phases is assumed.
r Low density and low viscosity of gas components result in fast gas dynamics that are on the order of 10 milliseconds. Since dynamics with time constants larger than 100 milliseconds are of interest in this study, quasi-steady state is assumed for gas phase transport.
r Water produced in the electrochemical reaction is assumed to be in vapor form.
r Gas mixture in both anode and cathode electrodes is assumed to be binary and Fick's law is used to describe gas diffusion in the diffusion media. This assumption is justified later in the paper.
r Since only the GDL is considered in the study, Knudsen diffusion is neglected.
r It is assumed that the velocity of the gas phase normal to the membrane is very low. Therefore, convective flow in the diffusion media and its effects on heat transfer are neglected. This assumption has been validated through scale analysis, where the Peclet number for the mass and heat transfer were found to be 10 −2 and 10 −3 , respectively. 46 Moreover, Benziger et al. showed that for lower oxygen molar fractions in the cathode (which is the case with air feed) the dominant mode of transport is diffusion. 47 Nevertheless, this assumption needs to be revisited based on the operating condition and cell configuration. For instance, Weber et al. found convective heat transfer to be significant. 23 r A small constant value is used for saturation at the GDL/GC interface to model the resistance to liquid flow at the surface. We acknowledge the need for a dynamic boundary condition based on the liquid saturation state inside the channel as was shown previously. 48,49 Nevertheless, this simple boundary condition provides a good starting point and was successfully employed in prior works. 19,35 r Liquid water in the channels is assumed to move with the gas phase as fine dispersed droplets (mist flow). Therefore a homogeneous model is used for pressure drop along the flow channels. Numerous studies have shown liquid water appearing in the channels in the form of a thin film or slugs. 50,51 We note that such effects in the channels can dramatically change the system's behavior. However, existence of liquid water in the channels is unfavorable for fuel cell operation and we assume operating conditions are chosen to avoid liquid water build-up in the channels.
With quasi-steady state assumption for gas phase transport, the dynamic states captured in the model include the temperature distribution, GDL saturation, and membrane water content. Next, each part of the model is presented separately. Flow channel model.-As mentioned earlier, a quasi-steady state is assumed for the reactant transport in the flow channels. The channel model is essentially a set of 12 ODEs that resolve 6 variables in each channel spatially in the along-the-channel direction. An additional ODE accounts for the variations in the coolant temperature along the flow channel. It is assumed that the contents within each element along the channel are well mixed and there is no gradient in the perpendicular direction.
A representative channel element is illustrated in Fig. 2. The nodal values of each element serve as boundary conditions for the sandwich model. Furthermore, the sandwich model determines the amount of water and reactant that is exchanged between the channel and GDL.
Hydrogen and oxygen flow rates in the anode and cathode gas channels are determined by: where the first term is due to reactant consumption (Faraday's law), the second term is due to reactant crossover through the membrane, and the third term represents the amount of O 2 (H 2 ) that reacts with the crossed H 2 (O 2 ) to produce water (see the membrane model section). A complete nomenclature and parameter list is provided in the Appendix. When a flow-through cell configuration is used in automotive applications, the anode gas stream is recirculated to ensure a high utilization of the hydrogen feed. In such cases, it is important to account for the nitrogen buildup in the anode gas stream since it dilutes the hydrogen in the channels. Therefore, the model includes nitrogen crossover in the membrane. Due to this crossover, nitrogen flow rates change along the flow channels: d N N 2 ,an (y) dy = wN N 2 ,cross (y) [ 3 ] d N N 2 ,ca (y) dy = −wN N 2 ,cross (y) [ 4 ] where N N 2 ,cross is the crossover nitrogen flux which is given by Eq. 32.
It should be noted that in order to account for the nitrogen crossover accurately, its partial pressures at both membrane interfaces are required. However, since nitrogen flow rates in the through-the-membrane direction are negligible, the partial pressures of N 2 in Eq. 32 are approximated by their respective values in the gas channels. The liquid water flow rate in the channels is determined by the amount of water vapor that condenses in the channels as well as the liquid water that penetrates the channels from the GDLs. This is described by: where the first term is due to condensation and the second term is due to the liquid water intrusion from the GDL. The saturation pressure is given by the following empirical fit: The water vapor flow rate is similarly affected by the phase change in the channel as well as the vapor exchange rate between the channel and the GDL and is resolved by: {P wv (y) − P sat (y)} + wN wv,GDL,e e ∈ {an, ca} [7] Note that in order to capture the phase change behavior in the channels correctly, a hybrid systems approach has to be utilized. Each channel has 2 modes corresponding to either single or two-phase flow. Overall, there are 4 modes for both channels depending on the flow conditions. To switch between these modes, functions similar to those in Ref. 52 are used. The remaining equations describe the temperature variations of each gas stream and the coolant: − T GC,an (y) C P,H 2 N H 2 (y) + C P,wv N wv,an (y) + C P,lw N lw,an (y) + C P,N 2 N N 2 ,an [8] dT GC,ca (y) dy C P,O 2 N O 2 (y) + C P,wv N wv,ca (y) + C P,lw N lw,ca (y) + C P,N 2 N N 2 ,ca [9] dT cool (y) dy = U cool A cool T 1 (y)+T 8 (y) 2 − T cool (y) C P,lw N cool [10] where indexed temperatures are evaluated at the interfaces depicted in Fig. 1.
To model the pressure drop along the channel, a homogeneous two-phase model is implemented. This is in line with the assumption of mist flow in the channels. In short, the pressure drop is given by: In the above relation, f is the friction factor and f Re is given by 53 : where α is the channel aspect ratio defined as: The two-phase mixture density is determined by: where x is the mass quality of the gas phase. Finally, the following model developed in Ref. 54 is used for the two-phase viscosity calculation, because it was shown to better match the experimental data 55 : In the above relations, gas mixture properties are obtained by molar averaging of individual gas constituent properties. Viscosity of individual gas components is calculated by Sutherland's formula 56 at the gas stream temperature along the channel.
To complete the pressure drop model, losses due to the entrance and exit regions are also added to the model. These regions correspond to the inactive parts of the cell area. The pressure drops in these regions are given by 57 : [17] where K c and K e are contraction and expansion loss coefficients that depend on the geometry, and K ∞ is the Hagenbach's factor, which, for rectangular channels, can be calculated by 57 : K ∞ = 0.6796 + 1.2197α + 3.3089α 2 − 9.592α 3 + 8.9089α 4 − 2.9949α 5 [18] Overall, Equations 1-11 describe the channel model in this study. Since the model is developed for a co-flow configuration, these equations can be solved together as an initial value problem with specified conditions at the inlet. If the model is desired to be solved for counterflow configuration, the problem will change into a boundary value problem, which can be solved with shooting methods, albeit at a higher computational cost due to the iterative nature of such algorithms.
Membrane model.-Traditionally, models developed for water transport across the membrane can be categorized as diffusivetype 58,59 and hydraulic-type 60 models. The former describes water transport based on the gradient in water concentration across the membrane and treats the membrane as a single homogeneous phase, whereas the latter relates transport to the difference in hydraulic pressure in the membrane, which is considered to be a two-phase system. Diffusive-type models are usually best capable of describing water transport in vapor equilibrated membranes. On the other hand, hydraulic-type models are better suited for fully humidified and liquid equilibrated membranes. A combination model based on concentrated solution theory was put forth by Weber and Newman, 61 which provides good explanations for the physical phenomena. Since computational efficiency is critical in our study, a simplified combination model is used, where water transport is affected by both the gradient in the water concentration and the difference in the liquid hydraulic pressure. Therefore, the membrane water crossover is given by: where positive direction is assumed to be from the anode towards the cathode. The first term is due to the electro-osmotic drag (EOD), and the second and third terms are back diffusion due to the concentration gradient and difference in the hydraulic pressure across the membrane, respectively. The last term is the water flux due to thermo-osmosis, which is shown to be from the cold to the hot side for hydrophilic membranes. 1 Notice that the last two terms are only present when the membrane is in contact with liquid water. Therefore, a switch (sw l ) is used that is turned on when the reduced saturation (see the gas diffusion layer model section) at one of the membrane interfaces becomes positive. Many empirical relations are reported for the EOD coefficient (n d ) and membrane diffusivity (D w,mb ) that vary greatly in some cases.
Similarly, values reported for the hydraulic permeability (K P,mb ) have large variations between 1.8 × 10 −18 m 260 and 2 × 10 −20 m 2 . 62 Here, a constant value of 1 is used for n d similar to Ref. 63.
The discrepancies between various membrane diffusivity values reported in the literature are significant and need careful treatment. Most models use the diffusivity values obtained by Zawodzinski et al. 64 Their diffusivity function goes through a maximum at a membrane water content of about 3, after which it shows a decreasing trend with the water content. It is worth noting that they first obtained the chemical diffusivity values and then used the Darken factor to transform them into the Fickian diffusivity values for the membrane. Many other experiments including those of Ge et al., 65 Ye et al. 66 and Motupally et al. 67 also show a peak in the diffusivity at water content values around 3. Most recently, Olesen et al. 68 showed that the peak in membrane diffusivity is purely a numerical artifact that is due to the change in the curvature of the membrane water uptake curve. The Darken factor, due to its logarithmic nature, is highly sensitive to this change in curvature. As such, lack of data around this inflection point could result in significant errors in diffusivity calculations. Using the chemical diffusivity of Zhao et al. 69 and the water uptake isotherms of Hinatsu et al., 70 Olesen et al. 68 found Fickian diffusivity values for Nafion membranes with no swelling that can be fit with the following sigmoid function 71 : [20] where λ trans = 2.6225 is the transition water content at which the diffusivity dependence on the water content switches from exponential to a second order polynomial. δ trans = 0.8728 is the transition interval width. If the membrane swelling is accounted for, we arrive at the following diffusivity relation 71 : D w,mb,no swelling [21] where D w,mb,no swelling is given by Eq. 20. These diffusivity values along with those of Zawodzinski et al. are shown in Fig. 3. It is interesting to note that Caulk et al., 72 through their permeation experiments and correctly accounting for all resistances to water transport, also observed similar trends in membrane diffusivity as those obtained by Eq. 21. This was also the case in experiments of van Bussel et al. 73 Moreover, this type of behavior is consistent with the observations of Zhao et al. that a first and second solvation shell form in the membrane as the  water content is increased. 69 Based on these arguments, we use Eq. 21 for membrane diffusivity in our model.
For the hydraulic permeability of the membrane, the following function is used that is obtained from data by Duan et al. 74 for Nafion N115: This polynomial fit is shown in Fig. 4. The data is not readily applicable to other membranes. However, in lack of other experimental data, this data is used as a starting point for the modeling purposes of this work. Similarly, experimental data of Kim and Mench 1 is used for the thermal diffusivity in a Nafion membrane. As noted earlier, there is a lot of uncertainty associated with these parameters in the literature. This is mainly due to the differences between the experimental setup and materials used in different tests. Therefore, the use of fitting parameters in Eq. 19 is essential to obtain a good agreement with experimental data. In fact, most of the studies that compare modeling results to experimental data take advantage of such parameters. 35,61 Therefore, selecting where to use fitting  parameters is as important as choosing the right empirical relation from the literature. We believe that a multiplier for the water back diffusion due to the concentration gradient is best suited to calibrate the presented model with experimental data.
Finally, the membrane water content dynamics is given by: which is a time-dependent nonlinear diffusion law. Traditionally, it was assumed that the water in ionomer phase inside the membrane is in equilibrium with the water outside of it. This assumption results in a Fickian diffusion of water across the membrane, where the only resistance to water transport is that due to diffusion. 58 However, recent experimental studies have revealed the existence of other resistances at the interface of the membrane equilibrated with water vapor, which are termed as interfacial resistances. 65,69,[75][76][77][78] Inclusion of such resistances has been shown to change water transport across the membrane dramatically. 75 It is important to note that these resistances are only present at the interface that is in contact with water vapor, because water uptake from a liquid reservoir has fast dynamics and quickly reaches the equilibrium state.
To model these interfacial effects, flux boundary conditions are required at each interface as given by: where k an and k ca are the mass transfer coefficients at the anode and cathode interfaces, respectively. These coefficients depend on both the water activity at the interface and the temperature. Kientiz et al. 79 provided an exponential relation for the mass transfer coefficient based on water activity alone. A temperature dependent relation was given by Ge et al., 65 which is used in this study: where f v is the volume fraction of water in the membrane and is given by: Eq. 26 is used for the absorbing interface and Eq. 27 is used for the desorbing interface of the membrane. When an interface is in contact with liquid water, the corresponding mass transfer coefficient is increased by three orders of magnitude to model the effects of vanishing interfacial resistance. Furthermore, λ eq an/ca in Eqs. 24 and 25 denotes the equilibrium membrane water content, while λ an/ca denotes the water content just inside the membrane. To find the equilibrium water content for a vapor equilibrated membrane, a linear interpolation was proposed by Ge et al. 65 based on the uptake data of Springer et al. 58 at 303 K and that of Hinatsu et al. 70 in contact with liquid water [30] To avoid the discontinuous jump when the membrane comes into contact with liquid water, a sigmoid function of the following form is used 80 : where h shape is a shape factor whose larger values result in a sharper transition between the two modes. A value of 8 is used in this model. Gas permeation through the membrane could also affect the cell performance. Most importantly, nitrogen buildup in the recirculating anode feed should be accounted for, which was addressed in the section on flow channel model. However, reactant crossover can also be important, since it results in performance loss. 81 The model accounts for the crossover of both reactants in addition to nitrogen, which was mentioned earlier. It is assumed that the hydrogen crossing over to the cathode reacts with the available oxygen immediately to produce water. The same assumption applies to the oxygen that crosses over to the anode side. Overall, the crossover fluxes are given by: where p j,ca/an denotes the partial pressure of gas species at the membrane interface and ψ j denotes the membrane gas permeability for the different species and is given in Table I. As mentioned earlier, the nitrogen partial pressure can be approximated by its respective value in the channel with relatively small errors. However, the partial pressure of oxygen and hydrogen should be evaluated at the membrane interface due to their non-negligible flux in the GDLs.
Reduced membrane model.-For state-of-the-art membranes that are thinner than 30 μm, the water content distribution is very close to a linear distribution. As such, the PDE in Eq. 23 can be replaced by an ODE in the form of 35 : ) unless CC License in place (see abstract This simplification results in significant computational advantages and is therefore used in this work. It should also be noted that compared to the original membrane model, only the time update (Eq. 23) is performed differently in the reduced model. A schematic of the simplified membrane model is shown in Fig. 5. Overall, there are four unknowns in the membrane model, namely, the equilibrium water contents at either side of the membrane (λ eq an/ca ) as well as the water contents inside the membrane (λ an/ca ). These unknowns are found by solving four equations that ensure the continuity of water flux across the different layers (see the overall water balance section).

Gas diffusion layer model.-
To complete the cell model, the twophase transport in the GDLs needs to be described. Gas phase transport is modeled by the steady-state Fickian diffusion: where R v is the molar phase change rate (given by Eq. 43) multiplied by the segment length and c j is the species concentration determined by: Note that nitrogen is neither consumed nor produced in any of the electrochemical reactions. As such, it is not considered in the diffusion model. The binary diffusivities are evaluated at the electrode pressure and average temperature using 9 : where D AB denotes binary diffusivity. Table II shows the diffusivities for each of the gas components with units of m 2 /s. These values need to be corrected to account for the porous nature of the GDL as well as liquid accumulation in the GDL, since the presence of liquid in the GDL effectively blocks the pores to gas phase transport. This correction is achieved through the following relation initially proposed by Nam and Kaviany, 86 but with the modified exponents based on Ref. 87: where s avg is the average liquid saturation level in the GDL and D 0 is the nominal diffusivity value calculated from Table II. Note that with this formulation, the concentration distributions of O 2 and H 2 are linear (steady-state Fick's law with constant diffusivity throughout the domain yields a linear distribution). As discussed earlier, multicomponent diffusion is required to describe the details of gas transport in the electrodes. However, in a pilot study we found that accounting for the multicomponent effects only becomes important at high current densities and such high values are usually avoided in applications due to the associated power losses. Kulikovsky 16 has also shown that using Fick's law introduces small errors while leading to much faster algorithms. Therefore, Fick's law is a good approximation for most applications.
The liquid phase movement is described by Darcy's law: where K rel = S 3 is the relative permeability and immobile saturation (also known as irreducible saturation) is used to define the reduced saturation 86 as follows: Finally, the capillary pressure is related to the saturation through the Leverette J-function: Using the surface contact angle of the GDL to describe its wettability and the application of the Leverette J-function that was originally developed for homogeneous sand is reported to be insufficient in describing the capillarity in a typical GDL material. 88 Other capillary relations are also reported in the literature. 89,90 Nevertheless, the above formulation is used in this study to enable a comparison of the results to those reported in the literature. The effects of incorporating other capillary pressure models is out of the scope of current work and should be investigated separately.
The two phases are coupled through the evaporation and condensation rates denoted by r v . A saturation concentration is calculated based on the temperature distribution in the GDL using Eqs. 6 and 36. The molar phase change rate is then given by 35 : where k evp/con is the volumetric phase change rate. It is suggested that a high value be used for this parameter to ensure that the twophase equilibrium assumption is not violated. 86 The selection of this parameter is important for a dynamic model, as it can affect the time scales of liquid accumulation and temperature distribution in the cell.
Here, a value of 100 is used in agreement with prior studies. 52 The above equation creates a nonlinear coupling between the two phases. To accurately account for this coupling, a fully coupled solution methodology must be employed, which is not amenable to realtime computations. To avoid this problem, we solve the gas phase transport without the nonlinear term due to phase change. This yields a spatially invariant water vapor flux as well as a linear distribution of the water vapor concentration in the GDL using Eq. 35. Using the water vapor distribution found this way, the source term is now calculated and used for the liquid phase transport (Eq. 40). Furthermore, the water vapor flux that is found initially is changed to account for the phase change: N wv,GDL,after phase change = N wv,GDL,before phase change + δ GDL 0 r v dx [44] This corrected flux value is used in Eq. 7.
Overall water balance.-With the simplified models for the membrane and the GDL described in the previous sections, the overall water balance can now be explained in detail. Based on the assumptions of this study for water transport (see the reduced membrane model and gas diffusion layer model sections), there is a constant (spatially invariant) water vapor flux on each side of the cell. The constant flux assumption results in very simple conditions for continuity across the layers.
At the GDL/GC interface, continuity implies that the water vapor flux at the interface is equal to that in the GDL: where h m = D h /(Sh D wv,e ) denotes the resistance to convective transport and is related to the channel hydraulic diameter (D h ), Sherwood number (Sh), and vapor diffusivity in free space (D wv,e ) given in Table II. A value of 2.693 is used for Sh, since the fuel cell channel is permeable on only one side. With this convective resistance, the water vapor flux between the membrane interface and channel can be directly described by: where R tot,GDL denotes the total resistance to water vapor transport between the membrane interface and the gas channel. This resistance consists of the convective resistance described above as well as the diffusive resistance in the GDL. That is: At the membrane interface, the following conditions are imposed for the water vapor flow rate: Using the above equation with Eq. 24 and Eq. 25, we obtain the water contents inside the membrane boundaries ( Fig. 5): [50] where k an/ca denotes the membrane mass transfer coefficient that was defined in Eqs. 26 and 27. Finally, from Eq. 19 we have: Using Eqs. 51 and 52 in Eq. 48 results in two algebraic equations that, along with the ODE given by Eq. 33, define a system of nonlinear differential algebraic equations (DAEs). The two algebraic unknowns of this system are the water vapor concentrations at the membrane boundaries, while the membrane water content (λ mb in Fig. 5) is the dynamic variable. The solution to this DAE will be discussed in the numerical implementation section.
Voltage model.-The cell voltage is determined by 91 : [53] where the first term is the open-circuit voltage (OCV), the second term is due to the cathode activation overpotential, third and fourth terms are due to the concentration overpotentials, and the last term describes the ohmic losses. The open-circuit voltage is given by: In the activation overpotential term, i cross accounts for the crossover current due to hydrogen permeation through the membrane which is given by 92 : To account for the voltage loss due to the catalyst coverage by liquid water, the cathodic exchange current density is determined by: where s cat denotes the saturation at the cathode catalyst layer and n v is a fitting parameter, 48 which is set to 1 in the present study, but can be used to calibrate the model with experimental data. The activation energy associated with oxygen reduction on platinum, E act , is 67 KJ/mol. 91 The reference pressure is 1 atm and the reference temperature is 303 K. Finally, the limiting current densities are found by setting the reactant concentrations at the catalyst layers equal to zero: in the denominator is due to the convective resistance at the GDL/GC interface. Note that free space diffusivity values should be used to model the convective resistance at the interface. Furthermore, RT mb H H 2 /O 2 is used to account for the reactant absorption into the ionomer phase, where H H 2 /O 2 denotes the Henry's Law constant for reactants. H H 2 is relatively insensitive to changes in temperature and a constant value of 4.5 × 10 −2 atm · m 3 mol −1 is used here. However, H O 2 is temperature dependent and is given by 20 : T mb + 14.1 [59] We note that use of Henry's Law in the above form results in a discontinuous jump in reactant concentrations at the catalyst site, which is attributed to species transport across a phase boundary. 20 Finally, the ohmic resistance consists of the electrical resistances of various layers, contact resistances, and membrane protonic resistance which is given by 58 : Eq. 53 is solved implicitly for a given cell voltage to obtain the current density at the desired location along the channel.
Temperature model.-The temperature dynamics are described by the heat equation: where the effective heat capacity and conductivities are found by volume averaging 1,91 :  [63] The following relation, suggested in Ref. 93, is used for the membrane conductivity: The source term is given by:  [65] where i ORR is the transfer current of the ORR and ORR denotes its Peltier coefficient 23 : T mb 298 [66] As shown in Eq. 65, only reversible and irreversible heat due to the ORR is considered, since the contribution of the HOR is usually negligible. It should be noted that since the catalyst layers are treated as interfaces in this model, their respective source terms are applied at the boundary nodes of the membrane. The source term for the membrane is applied on all membrane nodes. The heat of sorption is included in the model through the following relation 94 :  [67] Note that the heat released by water sorption is slightly more than the heat of condensation, because the water bound in Nafion is not in liquid phase. The same can be said for the heat absorbed through water desorption. Finally, H v is the phase change enthalpy, which is found by the following polynomial fit:

Numerical Implementation
In this section boundary conditions for various model equations are provided. Furthermore, numerical implementation of the model, which is critical for computational efficiency, is discussed.
Boundary conditions.-The channel model only requires the inlet values, which are provided as part of the operating conditions. The boundary conditions for the GDL saturation PDE, Eq. 40, include a zero flux at the membrane interface and a specified value of S = 0.001 at the channel interface. Boundary condition at the channel interface can be improved based on the liquid accumulation in the channel, 49 but this is out of the scope of the present study.
Flux relations similar to Eq. 45 are used for oxygen and hydrogen, where their fluxes at the membrane boundaries are set to their respective consumption rate. This determines the concentrations of all gas species throughout the GDL .
To find the temperature distribution in the cell, two boundary conditions are required. Flux boundary conditions are used on both sides, where it is assumed that the cell heat is dissipated by convection with the gas streams in the channels: [69] Similarly, convective heat transfer with coolant and gas flow in the channels is taken into account for the bipolar plates. It is noted that Bhaiya et al. 6 used Dirichlet boundary conditions for the cell temperature at the GDL/GC interface. However, the convective heat transfer is important in this case and should be included in the temperature model. We believe a mixed boundary condition is best suited for this interface to account for the heat dissipation through both the gas channel and the land. Alternatively, one could use the conformal mapping approach presented by Weber et al. 43 to capture the effects of landings. We use only the flux boundary condition in our study, since the lands and their effects are not considered in the model. Fig. 6 illustrates the pseudo-2D framework as well as the numerical grid used in the through-the-membrane direction. Fig. 7 shows the specific solution method employed in this work. The distinctive feature of the approach is that dynamic sub-models are solved separately, which significantly reduces the computational requirements. At the heart of this approach lies the assumption that the main coupling between such dynamics can be captured through boundary conditions and required source terms. Therefore, fully coupled solutions are not necessary. The advantages of this approach include:

Solution method.-
r Sub-problems can be solved using different numerical schemes that are best suited for the specific problem.
r Parallel computing can be utilized, since time updates are performed separately.
In particular, the backward Euler method (fully implicit) is used for the temperature PDE (Eq. 61), while the saturation PDE (Eq. 40) is solved using the Beam and Warming implicit method. 95 It should be noted that micrometer scale of the fuel cell layers results in numerical stability issues for explicit solution techniques. Therefore, with such explicit methods, very small time step sizes are required to ensure that the stability criteria are satisfied. These issues are avoided in this work through implicit methods, which in turn enable the use of larger time step sizes and reduce the computational time needed to simulate the model. The cell sandwich model solves for both the current density and water vapor distribution in the cell. It also updates the membrane water content at each node. To find the current density at each nodal point along the channel, Eq. 53 is solved for a specified cell voltage using the inverse quadratic interpolation protected by bisection to ensure the robustness of the solution. 96 The water distribution is found by solving a DAE system as described in the overall water balance section. This DAE system is solved using Newton's method, where the water vapor concentration in the channels and the previous membrane water content provide a good initial point for the solver. Furthermore, since the unknown variables differ by several orders of magnitude, the Jacobian matrix has a very high condition number. This issue can be addressed through scaling the variables by their nominal values.
The computational grid consists of 25 nodes in the along-thechannel direction, 20 nodes for each GDL and 3 nodes for the membrane in the through-the-membrane direction. This grid size was chosen to minimize grid dependency while maintaining computational efficiency. Generally, fewer nodes in the along-the-channel direction reduces computational time while the grid size in the through-themembrane direction does not have a significant effect on computational efficiency. However, the number of nodes along the channels has to be high enough to ensure convergence. In particular, numerical stiffness of the channel model (caused by simultaneous solution of the temperature and gas concentration distribution) limits the number of nodes in the along-the-channel direction to above 20. This limit in the  number of nodes may change with operating conditions, meaning that for some conditions convergence can be achieved with fewer nodes along the channel. In this study, 25 nodes were used along the channel, since it was found to result in convergent solutions in all the test cases. As for the time step size, a time step of 100 milliseconds is used in this work. Smaller time step sizes should not be used with this model since gas dynamics (that are neglected in the current model) would be present at those time scales. On the other hand, larger time steps may be used, provided that the solution does not show oscillations. In our simulations, we observed stable and convergent solutions with time steps up to 200 milliseconds. Nevertheless, we chose 100 milliseconds to obtain better accuracy.
In terms of computational time, the 800 second simulation results provided in the case study on transient behavior below were obtained in 431 seconds with a MATLAB implementation of the model on a personal computer (2.4 GHz CPU). Moreover, the dynamic cycles presented in the voltage and relative humidity cycling section that simulate 1200 seconds were obtained in an average time of 611 seconds (with the longest computation time of 642 seconds). Therefore, the model is almost twice faster than real time. Furthermore, an initial compiled version of the code shows a computational performance improvement of up to four times. To the best of our knowledge, this is the first time that a model with such computational efficiency incorporating this level of physical details is reported in the literature.

Results and Discussion
Transient behavior -case study.-The transient behavior of the fuel cell is investigated in this section through a case study. The operating conditions are described in Table III. It is assumed that the cell operates under a potentio-static mode; i.e., a cell voltage is given as an input to the model and current density is one of the model outputs.  Table III: a) average current density b) average membrane water content c) average cell temperature d) average GDL saturation.

Initial conditions
Average system dynamics are illustrated in Fig. 8. It is evident from the figure that the slowest transient is due to liquid water accumulation in the GDLs, while the membrane water uptake and temperature both have faster dynamics with time constants on the order of 10 seconds. The cell current density is affected by all of these transient behaviors.
In particular, during the initial 30 seconds start-up, both the membrane water content and cell temperature increase rapidly. Better hydration of the membrane lowers the ohmic losses in the cell, thereby increasing the current density. On the other hand, the effects of the temperature increase are twofold; first, it results in a decline in the thermodynamic cell voltage, which in turn tends to lower the current density. Second, elevated temperatures improve the reaction kinetics (Eq. 56), thereby improving cell performance. The temperature also has indirect effects on the cell performance by affecting the membrane conductivity; temperature increase improves the membrane ionic conductivity (Eq. 60) at a given water content level, but generally tends to lower the relative humidity of the gas phase adjacent to the membrane, which in turn affects the membrane water uptake. Under normal operating conditions, temperature increase usually results in a higher ionic resistance of the membrane. For the presented results, it is seen that during the initial start-up. better membrane hydration has resulted in an increase in the current density.
The sudden increase in the current density is then followed by a slow decay due to the GDL and cathode catalyst layer flooding. It is observed that this decay occurs over a longer time period compared to the initial transients. Furthermore, the figure shows a shoulder in the current density dynamics, where the slope of the decay is further decreased at about 250 seconds. This corresponds to the time that the cathode GDL saturation has reached close to 90 percent of its steady state value. At this point the rate of flooding is decreased due to the lowered current density and water production. At about the same time, liquid water starts to appear in the anode GDL. However, anode flooding does not restrict the system performance considerably, since hydrogen diffusion in the anode is fast. Therefore, the decay in the current density continues at a lower rate after 250 seconds. As for the cell temperature, it closely follows the current density dynamics, since the main heat source is due to the reversible and irreversible heat generation from the ORR, which is directly related to the current density.
The current density distribution and channel water flow rates are shown in Fig. 9. In the along-the-channel direction, the current density first increases slightly due to the membrane hydration. At about 50 percent of the channel length, the current density drops significantly due to the GDL and catalyst layer flooding. After this point, the current density decreases further towards the channel outlet, which is attributed to reactant dilution. As for the water flow rates, the vapor flow rate in the cathode increases initially due to water generation and the positive net water transport from anode to cathode. However, towards the channel outlet the flow rate remains relatively constant due to lower water generation as well as higher back diffusion to the anode. The higher back diffusion close to the channel outlet can also be observed from the increase in the anode vapor flow rate. Finally, it is seen that liquid water starts to appear in the channels at about 0.6 and 0.8 of the distance from the cathode and anode inlets, respectively. It should also be noted that the flow rate of liquid water is one order of magnitude smaller than that of water vapor. This is in part due to the homogeneous flow assumption that is used in the channel model.
Finally, Fig. 10 shows the temperature distribution and GDL saturation in the cell at time instants of 50, 200, and 800 seconds. Note that initially, the temperature increases monotonically from the channel inlet to its outlet. This is due to the assumption that the coolant flows in the same direction as the reactants. Therefore, towards the outlet the coolant temperature is higher and not effective for heat removal. As time goes on, the location of the maximum temperature moves towards the middle of the channel. This is attributed to the flooding in the cathode GDL, which results in lower current densities towards the channel outlet. Therefore, after the flooding occurs in the GDL, the temperature increases from the inlet to the onset of the GDL flooding. At the point of the GDL flooding, the temperature drops due to a lower current density. After this point, the temperature increases slightly towards the channel outlet due to the increasing coolant temperatures. Another important observation is the direction of the liquid water flow in the GDLs. In the cathode GDL liquid water flows from the membrane towards the channel. However, in the anode GDL this direction is reversed, meaning that liquid water moves towards the membrane. This is known as the phase-change-induced (PCI) flow. 23,27,97 This type of flow is due to the temperature gradients in the cell that cause a gradient in the vapor saturation pressure. Therefore, water vapor moves along the temperature gradient and condenses on the cooler side of the cell. The condensed water then moves back towards the membrane by capillary flow, creating a recirculation mechanism for water on the anode side. Also note that this does not usually happen on the cathode side in spite of temperature gradients, which is due to the slower diffusion of water vapor in the cathode GDL as well as water production on the cathode side. An important implication of this phenomena is that the anode side of the membrane remains dry in spite of the liquid water presence in the anode GDL. This in turn results in increased back transport of water from the cathode to the anode side. To the authors' best knowledge, it is the first time that such effects are captured in a computationally efficient model.
Water transport in the cell.-Different water transport mechanisms are investigated in this section to determine their contribution to the overall membrane water crossover. To this end, each of the terms in Eq. 19 is plotted in Fig. 11. The values are obtained by averaging in the along-the-channel direction. In the figure, positive sign denotes water flux from the anode to the cathode side. The EOD flux follows the current density dynamics. The diffusive transport is negative on both sides, where the initial transients correspond to the membrane water uptake. These two fluxes settle to a negative steady state value, meaning that the net diffusive water flux is towards the anode as expected. The thermo-osmosis flux and hydraulic fluxes are initially zero, since no liquid water exists on the membrane interfaces. As liquid water builds up in the GDLs and the reduced saturation becomes nonzero, both terms start to appear, with the thermo-osmosis flux being from the anode to the cathode, whereas the hydraulic flux is in the opposite direction. However, note that the water fluxes due to these terms are orders of magnitude smaller than the other fluxes. Therefore, the water flux across the membrane is dominated by the competing effects of the EOD and diffusive transport; this is explained further below. In this case study, the diffusive transport is more significant, resulting in a net water transport from the cathode to the anode side. We note that this is in agreement with experimental results, where significant water removal from the anode side of the cell was observed. 98 As pointed out by Duan et al., 74 experimental data suggests that the water flux due to the hydraulic pressure is negligible compared to that due to the diffusion by concentration difference. Therefore, when a concentration difference is present, the water flux by hydraulic pressure can be neglected. This is often the case for membrane humidifiers. However, in an operating fuel cell, the membrane is usually in contact with liquid water on both sides, since both gas feeds are humidified. This results in similar water activity on both of the membrane interfaces, thereby reducing the diffusive flux significantly. Duan et al. argue that in such cases, the hydraulic transport can become important compared to the diffusive flux. Nevertheless, we note first that the diffusive flux in our model is based on the difference in the membrane water content, which can be present even if both electrodes are flooded. This is a direct result of the combination model, where the equilibrium water content is found by interpolation between the fully vapor saturated and fully liquid saturated cases. In particular, the saturation state of the cathode GDL is usually higher than that of the anode GDL. Therefore, even if both electrodes are flooded, the diffusive flux can still be significant. Secondly, in the special case that the diffusive flux is indeed insignificant, the flux due to the EOD is still an order of magnitude higher than the hydraulic flux. Therefore, we believe that within our modeling approach both hydraulic and thermoosmosis fluxes can be safely neglected for a large range of operating conditions. Lastly, we note that even though thermo-osmosis does not contribute significantly to the water balance in the cell, capturing the temperature gradients is absolutely necessary for a correct prediction of the water state in different cell layers, because such gradients can result in redistribution of water (e.g., PCI flow).
To have a clear picture, water crossover through the anode interface of the membrane is plotted against normalized channel location for three time instances in Fig. 12. It is observed that close to the channel inlet, water flows from the anode to the cathode due to the EOD. At about 15 percent of the channel length, this direction is reversed as the diffusive back flux becomes significant enough to overcome the EOD. It can also be seen that the overall flux close to the channel inlet does not change from 50 to 800 seconds, since there is no flooding in that region and transient behavior is fast. However, the flux close to the channel outlet is affected by liquid accumulation, which in turn affects the membrane water content distribution. At the location  where the cathode GDL floods, water back transport is enhanced. As mentioned before, this is a direct result of the PCI flow, where liquid water in the cathode GDL forms on the membrane interface, whereas on the anode side liquid appears close to the channel first, leaving the membrane interface dry. As more water is transported to the anode, the water vapor concentration increases on the anode side, which in turn decreases the back transport. This explains the local minimum in the membrane water crossover along the channel. Fig. 13a shows the fraction of product water removed via the anode during cell operation with the same operating conditions as given in Table III. This fraction is defined as: where a positive value indicates water absorption into the membrane at the anode interface. This fraction is an important parameter in fuel cell operation, since it determines the hydration requirements for the anode feed. It is observed that the fraction is initially positive due to the membrane hydration. During this period, the membrane absorbs water from both the anode and the cathode. However, after the initial transient response the anode side of the membrane becomes a desorbing face on average, which makes β negative. As time goes on and the cathode GDL becomes saturated with liquid water, a liquid vapor permeation (LVP) mechanism 99 emerges that enhances water back transport, resulting in lower values of β. Fig. 13b illustrates the steady state values of β for different inlet relative humidities (relative humidity is the same for the anode and cathode streams). It is seen that increasing the inlet gas relative humidity results in more water to be transported to the anode (therefore a decrease in β). This can be explained by the better hydration of the membrane, which results in higher diffusivity values for water. Therefore, the water back transport towards the anode is enhanced with a better membrane hydration (also note that the EOD coefficient is kept at a constant value of 1 in our model and does not change with the membrane hydration). It is also worth pointing out that the values of β obtained here are on the higher end of the values reported in the literature, which is attributed to the small thickness of Nafion membrane used in the simulations (we have used 30 μm as the membrane thickness in our simulations while values reported in the modeling literature are usually greater than 50 μm). Furthermore, Cheah et al. 75 showed that the effective EOD coefficient can be an order of magnitude smaller than the nominal value due to the membrane interfacial resistances. Since those resistances are included in the model, this reduction in the effective EOD is captured and a higher water flux towards the anode is observed. Finally, the steady state water production and its removal rates via the anode and the cathode are shown in Fig. 14, which indicates that the amount of water removed through the anode initially increases with the current density before decreasing at higher current densities. The peak occurs at a current density of about 0.6 A/cm 2 . The decrease at higher current densities is attributed to a more significant EOD. This result agrees qualitatively with the experimental results of Thomas et al. 100 Hickner et al. 25 also observed similar behavior in their study, where the amount of liquid water in the anode was found to initially increase with the current density before decreasing at higher current densities.
Effect of stoichiometric ratio.-To investigate the effects of the stoichiometric ratio on the liquid water build-up in the cell, a sensitivity analysis is carried out, where ζ H 2 (ζ O 2 ) is varied while ζ O 2 (ζ H 2 ) is kept constant; other operating conditions are the same as Table III. The locations where the channel gas feeds become saturated at the steady state are shown in Fig. 15 for different stoichiometric ratios. Note that these locations are closely related to the location of liquid accumulation in the GDLs. It can be seen that an increase in ζ generally results in unsaturated gas feeds for a longer portion of the channels. In particular, the saturation of the anode feed is highly sensitive to the stoichiometric ratio of both sides, whereas that of the cathode feed is insensitive to the hydrogen stoichiometry. This is due to the fact that water back-flow towards the anode is highly dependent on water concentration in the cathode channel. As the cathode feed gets saturated further down the channel, the anode feed saturation is also delayed. This is due to the fact that water vapor in the anode channel comes from either the inlet or back transport through the membrane. As such, the anode channel saturation is highly dependent on the cathode saturation through the membrane water crossover. Due to water generation on the cathode side, the saturation in the cathode channel is not sensitive to the anode saturation. To have a better understanding of this behavior, β values are also plotted on the same figure. It is seen that β increases slightly with the anode stoichiometry, while it decreases significantly with cathode stoichiometry. This can be explained by the ability of each gas stream to hold water in vapor form, which increases with the gas stream stoichiometry (and therefore the gas flow rate). Therefore, increasing the stoichiometry ratio for the anode (cathode) gas stream results in a higher water removal rate through the anode (cathode) side of the cell. This is explained in more detail in the section on circuit analogy for water vapor transport below.
It is also observed that the anode gas feed becomes saturated further along the channel compared to the cathode feed. This is due to the fact that most of the water produced in the cell is carried out by the cathode gas stream. However, note that the molar fraction of hydrogen is about three times higher than that of oxygen in the cathode channel. As such, hydrogen depletion in the anode channel has a more significant impact on the saturation of the anode stream. This in turn mitigates the first effect, resulting in the fact that both gas feeds become saturated at about the same location along the channel as long as ζ H 2 is not considerably higher than ζ O 2 .
Voltage and relative humidity cycling.-In this section, cell behavior under dynamic cycles is investigated. Such dynamic cycles are usually encountered in automotive applications and the model has to be capable of performing these challenging simulations in a computationally efficient manner to be of practical use. Cell performance under voltage cycles with different inlet relative humidities is shown in Fig. 16. In all cases the high voltage (HV) is 0.90 V and the low voltage (LV) is 0.75 V. The sequence for each cycle is HV-LV-HV-LV and each cell voltage is maintained for a period of 300 seconds. All other conditions are the same as in Table III. It is observed that an increase in the inlet RH from 45% to 65% improves performance due to better membrane hydration. However, further increase in RH from 65% to 85% results in loss of performance due to the flooding in the cathode GDL. It is worth noting that even with 65% inlet RH there is some level of flooding in the cathode GDL, but the rate of flooding is very slow and its effect is not easily observed in the 300 second window. It is also seen that an increase in the current density is followed by a better membrane hydration. This means that under these conditions, higher water production rates (which results in a better membrane hydration) have a more significant influence on the membrane water content compared to the EOD (which tends to dehydrate the anode side at higher current densities). Fig. 17 illustrates the average cell dynamics under various relative humidity cycles. In particular, three simulations are performed with the RH of both the anode and cathode streams subject to dynamic cycles between low RH (LRH) and high RH (HRH) values: (1) HRH of 60% and LRH of 45%, (2) HRH of 75% and LRH of 45%, and (3) HRH of 90% and LRH of 45%. In all simulations the cycles follow HRH-LRH-HRH-LRH sequence and each condition is maintained for 300 seconds, allowing enough time to observe all the important dynamics. The cell voltage is maintained at 0.75 V for all cycles and all other conditions are the same as in Table III.
In case (1) the gas streams are only moderately humidified at HRH. Therefore, the increase in the RH results in better membrane humidification and lowers ionomer resistance, while no flooding is observed. As a result, the average current density is higher at HRH and lower at LRH. As for case (2), it is observed that an increase in the RH results in an overshoot in the current density; an initial increase in the current density due to better membrane humidification followed by a slow decay attributed to the GDL flooding as can be seen in Fig. 17d. When the RH is decreased, an undershoot in the current density is evident. This undershoot can be explained by the residual liquid water in the GDLs and the associated time constant is related to the liquid water evaporation rate in the GDL. Case (3) also shows the same behavior. However, the overshoot and undershoot are more pronounced due to the higher humidities at HRH in this case. Interestingly, at the end of HRH in case (3), the flooding effects completely overcome the humidification benefits and the current density is seen to be lower than its value at LRH.
In all the cases increasing the RH results in higher membrane water contents and more GDL flooding as expected. Note that only the cathode GDL saturation is shown, since saturation in the anode GDL was negligible in all the cases. Finally, temperature dynamics in all the simulations presented in this section follow the current density dynamics as was explained in the case study on transient behavior.
Circuit analogy for water vapor transport.-Based on the various transport resistances for water vapor movement in the cell, an equivalent electrical circuit model is constructed in this section to help with clarifying the relative importance of such resistances throughout the cell operation. The circuit is illustrated in Fig. 18. In this circuit, the water flow rate is given by the current, whereas the water concentration is defined as the voltage. Note that on each side of the membrane there exists three resistances between the membrane and   the gas channel; namely, a membrane interfacial resistance denoted by R mb,int , a diffusive resistance in the GDL denoted by R GDL , and a convective resistance at the gas channel interface denoted by R GC . The membrane also adds an additional diffusive resistance between the anode and cathode sides of the cell; this resistance is denoted by R mb,diff . Finally, the hydration state of the membrane is captured by a capacitance (C mb ), and a current sink on the anode side and a current source on the cathode side are used to model the EOD and water production, respectively. The water is produced on the cathode side of the membrane. As such, diffusing through the cathode GDL and exiting via the cathode gas stream is a more favorable flow path for water vapor, since this path has a lower resistance. In particular, there are three additional resistances to water back-flow, which include the membrane interfacial and diffusive resistances. Therefore, during the initial start-up period, most of the product water is removed via the cathode. However, as the cathode channel becomes saturated, the potential difference between the cathode catalyst layer and cathode gas channel vanishes. Therefore, water vapor can only diffuse through the membrane and flow towards the anode in spite of the higher resistance in this path. Additionally, water vapor that is now trapped in the cathode GDL condenses and forms liquid water on the membrane interface. With this liquid water formation, the membrane interfacial resistance vanishes, which reduces the overall resistance to water back-flow. In short, as the operation goes on, two mechanisms help enhance the water back-flow; namely, vanishing of the driving potential on the cathode side and reduction of the resistance to flow towards the anode side. This also explains the coupling between the location of cathode and anode gas feed saturation that was shown in Fig. 15. Specifically, as water vapor concentration increases in either of the channels or they become saturated, the potential difference driving water vapor towards that channel decreases. This in turn results in more water flow towards the unsaturated channel. Effects of reactant depletion can be described similarly, since it essentially translates into higher concentration of water vapor. Therefore, we believe that such equivalent circuit models, through their simplistic nature, can be of great utility in illustrating the overall water balance in the cell.

Conclusions
This study is aimed at filling the existing gap between high fidelity and low computational load in modeling polymer electrolyte fuel cells. Models that can provide high fidelity in a computationally efficient way can play a pivotal role for online health monitoring purposes, where real-time simulations provide important information about system states that can be used to avoid critical operating conditions.
To this end, a two-phase non-isothermal cell-level model is presented for polymer electrolyte fuel cells with straight flow channels. The model accounts for both along-the-channel and through-themembrane variations. Noting the disparate time scales of the fuel cell dynamics as well as the high aspect ratio of various layers, spatiotemporal decoupling is utilized to formulate the model in a compu-tationally efficient manner while maintaining fidelity. In particular, such decoupling enables the use of efficient numerical schemes that are developed specifically for each individual sub-problem, thereby contributing to the overall efficiency of the model. The current implementation is about twice faster than real time in an interpreted computer language, and further improvements can be achieved by using more efficient compiled languages that are generally used in embedded applications.
The model is used to perform a case study on overall system dynamics, where different performance loss mechanisms are investigated. It is found that for humidified anode and cathode gas feeds, which is usually used in automotive applications, cathode flooding results in the most significant performance loss with a relatively slow and two-stage time constant. Furthermore, water transport mechanisms are investigated to determine contribution due to each individual mechanism to the overall membrane water crossover. It is observed that the electro-osmotic drag and diffusive transport are the dominant contributors to the overall water balance for typical operating conditions, whereas hydraulic permeation and thermo-osmosis are found to be insignificant. Moreover, it is found that for state-of-the-art membranes, a higher fraction of product water is removed via anode due to enhanced back transport.
The model is shown to successfully capture some of the recently observed phenomena including temperature driven flow. It is found that, in accordance with published experimental results, liquid water distribution throughout the cell is determined by the interplay between local current density, local heating and reactant depletion along the flow channels. The model is also shown to efficiently simulate dynamic cycles that are encountered in automotive applications. Finally, an electrical circuit analogy is provided to elucidate various resistances to water transport. It is suggested that the cathode channel saturation forces the water to flow back towards the anode in spite of the higher flow resistance in that path.
To summarize, the contribution of this work is threefold. First, it provides a comprehensive non-isothermal and two-phase transient model for PEM fuel cells that incorporates most of the phenomena that have been previously shown to affect the overall fuel cell performance. Second, through a spatio-temporal decoupling it provides a computationally efficient solution that can be used for online applications. To the best of the authors' knowledge, this is the first model with this level of fidelity and computational efficiency that allows for realtime simulations. Finally, water transport across various cell layers is elucidated based on the modeling results and relative importance of each transport phenomena is discussed.
Ongoing efforts in our group are focused at investigation of microporous layer effects on the overall water balance in the cell, probing liquid water accumulation inside the gas channels, as well as incorporating the land effects through the use of effective material properties and interfacial resistances. Phase change rate in gas channels 100 s −1 U g Gas convection coefficient 2 × 10 3 W · m −2 · K −1 U cool Coolant convection coefficient 6 × 10 3 W · m −2 · K −1 C P,O 2 Specific heat of oxygen 925 J · kg −1 · K −1 C P,H 2 Specific heat of hydrogen 14400 J · kg −1 · K −1 C P,N 2 Specific heat of nitrogen 1041 J · kg −1 · K −1 C P,lw Specific heat of liquid water 4180 J · kg −1 · K −1 C P,wv Specific heat of water vapor 1840 J · kg −1 · K −1 ρ Density -kg · m −3 μ Viscosity -kg · m −1 .s −1 x Mass quality of gas phase in mixture