A Mathematical Model toward Real-Time Monitoring of Automotive PEM Fuel Cells

A computationally ef ﬁ cient model toward real-time monitoring of automotive polymer electrolyte membrane (PEM) fuel cell stacks is developed. Computational ef ﬁ ciency is achieved by spatio-temporal decoupling of the problem, developing a new reduced-order model for water balance across the membrane electrode assembly (MEA), and de ﬁ ning a new variable for cathode catalyst utilization that captures the trade-off between proton and mass transport limitations without additional computational cost. Together, these considerations result in the model calculations to be carried out more than an order of magnitude faster than real time. Moreover, a new iterative scheme allows for simulation of counter- ﬂ ow operation and makes the model ﬂ exible for different ﬂ ow con ﬁ gurations. The proposed model is validated with a wide range of experimental performance measurements from two different fuel cells. Finally, simulation case studies are presented to demonstrate the prediction capabilities of the model.

Cost reduction and durability enhancement remain at the forefront of research and development efforts aimed at commercializing PEM fuel cell vehicles. Remarkable improvements have been achieved on both fronts over the past decades by significantly reducing the amount of precious metal catalysts used, while enhancing the durability of various stack components. 1 Thus far, efforts have been mostly focused on materials research to develop new material solutions or improve the existing ones in terms of performance and durability. Nevertheless, innovative system design and control solutions also play a significant role and hold a great promise in addressing these challenges, as has been demonstrated by the commercialized Toyota Mirai. 2 Therefore, novel fuel cell management systems must be developed alongside the new material developments in order to achieve widespread commercialization of fuel cell vehicles. Such management systems may also prove helpful in utilizing the full capability of current materials.
Similar to the state-of-the-art battery management systems that rely on reduced-order models to make critical decisions that affect the performance and lifetime of the battery pack, 3 a fuel cell management system depends on a model of the fuel cell stack and the auxiliary components. 4,5 This is in part necessitated by the fact that direct measurement of many of the critical variables on-board a vehicle is difficult, costly, and in many cases impossible. This motivates the need for real-time models that can run efficiently onboard the vehicle and act as software sensors, providing information about the distribution of critical variables in the fuel cell stack.
Since the seminal work of Springer et al., 6 fuel cell modeling has attracted significant attention, resulting in numerous models that vary in terms of the dimensionality, incorporated physics, and computational costs. In particular, earlier models were mostly lumped parameter with no distribution, 4 or 1D, capturing the distributions in either the through-the-membrane 7 or along-the-channel 8 directions. More recently, most of the models have evolved to 2D computational domains, capturing land-channel variations [9][10][11] or along the channel distributions. [12][13][14] Pseudo 3D 15,16 and 3D 17,18 models have also become popular. Additionally, the transient phenomena have been incorporated in many recent models. 9,11,12,14,19 In most cases, the evolution to higher dimensions and incorporation of transient behavior are accompanied by a higher complexity of the modeled physics. More specifically, accounting for multi-step oxygen reduction reaction (ORR) kinetics 20 and intricate two-phase flow and nonisothermal phenomena [9][10][11]18 has become common in recent literature. This ever-increasing complexity of the models has shed some light into the transport phenomena in the porous layers of the cell that were previously unidentified. 21 Moreover, it has enabled a closer look at some of the degradation issues that limit the lifetime of the stacks. 22,23 Nevertheless, this all comes at a significant computational cost that inhibits the use of these models for on-line estimation and diagnostics. On the other hand, the available computationally efficient models typically disregard or oversimplify the underlying physics to the point that they miss many critically important phenomena, such as the twophase flow 24 and non-isothermal effects. 25,26 Despite recent efforts, there remains a gap in the fuel cell modeling literature for models that strike a balance between fidelity and computational cost.
To fill the above-mentioned gap, the objective of this work is to develop a PEM fuel cell model that is geared toward real-time simulation on-board the vehicle and provides critical information for fuel cell management and control systems. To this end, we utilize the disparate time scales of the involved transient phenomena as well as the large aspect ratio of the cell components to spatio-temporally decouple the problem and obtain a computationally efficient solution. Development of a new reduced-order model for water balance across the cell and a simplified method to account for mass and charge transport trade-offs in the cathode catalyst layer further contributes to the computational efficiency of the model.
Here, our specific focus is on capturing temperature and water distributions throughout the cell, which is achieved through a pseudo-2D bi-domain modeling approach. In particular, this work draws from our earlier effort in this area 13,27 and enhances it as follows: i) the bi-domain modeling approach allows for capturing of the in-plane distributions, ii) the microporous layer (MPL) is explicitly accounted for and is no longer lumped with the gas diffusion layer (GDL), while the catalyst layer (CL) is treated as a component with finite thickness rather than an interface, iii) the counter-flow configuration is modeled while maintaining real-time computation capabilities, iv) the model is more efficiently implemented to allow for significant savings in computation time, and v) operating conditions. These modifications render the model suitable for real-time monitoring of unmeasured and critical states within the fuel cell stack.
The rest of the paper is organized as follows. The model formulation is presented first, describing the steps for developing the reduced-order model. The presented model is then extensively verified with experimental data from two stacks and shown to provide good agreements with the data with reasonable parameter values. Finally, some simulation case studies are presented to further demonstrate the model's capabilities, followed by a brief summary and concluding remarks. A complete list of the model parameters and variables and their corresponding descriptions is provided in the Appendix Table A·I.

Model Development
Modeling domain.-The modeling domain is shown in Fig. 1. This is a pseudo-2D bi-domain (P2D 2 ) model, which captures the distributions in all three physical directions in a decoupled manner. In particular, a 1D through-the-membrane transport model is developed for the MEA. This 1D model is solved twice, once under the channel and once under the land area.
The bi-domain modeling approach was first proposed by Zaglio et al. 28 to efficiently capture the in-plane distributions without the need for a 2D model. Simplifying the 2D modeling domain to a 1D domain through conformal mapping has also been proposed in the literature. [29][30][31] While this method benefits from the highest computational efficiency, it would not allow one to distinguish between the local conditions in different parts of the cell, which can be critical, especially within the context of degradation studies. Therefore, we use the bi-domain approach here to balance the computational cost and model fidelity requirements.
The underlying idea is that an effective transport length for each of the transport equations can be defined on either domain, i.e. under the channel or under the land area. This effective length can be used along with appropriate boundary conditions to solve the corresponding transport problems for heat, mass, and charged species. This is illustrated in Fig. 1. In particular, the effective mass transport path is longer under the land area, whereas the regions under the channels have longer charge and heat transport paths. Here we account for the differences in mass and heat transport paths, while the charge transport problems are treated similarly in both regions, since the corresponding differences are negligible. In each case, the effective transport lengths between the channel and the CL are calculated by taking the in-plane (IP) and through-plane (TP) transport lengths into account (see Fig. 1b where d TP/IP denotes the TP/IP transport length, δ DM is the sum of the thicknesses of the diffusive layers (CL, MPL, and GDL), w CH/LN is the channel/land width. The effective length defined above also accounts for the material anisotropy through the k TP/IP factor, which is the ratio of the relevant effective transport property in the TP direction to that in the IP direction. This factor is determined based on the anisotropic material properties for each transport problem, which include effective thermal conductivity, diffusivity, and permeability for the heat, gas, and liquid transport problems, respectively. Finally, the transport length at any intermediate location within the diffusive layers can be found by scaling the above equation accordingly. The bi-domain MEA model is then extended along the flow channel using a pseudo-2D (P2D) modeling approach, where the only coupling between the channel and MEA models are obtained through the channel boundary conditions that account for mass and heat exchange between the GDLs and flow channels (see Fig. 1c). 13 This approach is suitable for PEM fuel cell modeling due to the length scales in the along-the-channel and through-the-membrane directions.
Through-the-membrane transport model.-Here we describe the through-plane transport model. We use common notation in our formulation and the complete description for the model variables and parameters is provided in the Appendix Table A·I.
Ionomer water uptake and transport.-The dynamics of water content in the ionomer (λ) are governed by different transport mechanisms as well as chemical (due to direct hydrogen combustion) and electrochemical (due to ORR) water generation. In particular, the water dynamics in the ionomer phase are given by three ordinary differential equations (ODEs): Note that thermo-osmosis is shown to drive water from the cold to the hot side for a hydrophilic membrane. 32 As a convention, a positive flux denotes water flow toward the cathode. The membrane water transport properties are given in Table I. There is considerable uncertainty surrounding such transport properties, 33,34 which warrants use of fitting parameters for membrane water flux, when agreement with experimental data is desired. For this reason, we have introduced a scaling factor for water diffusion, ξ diff,mb , which is identified during model parameterization. We also note that a similar model tuning strategy has been suggested in the literature. 35 In Eqs. 2-3, the water absorption/desorption source term, S ad , is given by: where λ* is the dynamic quasi-equilibrium water content for the ionomer used to model the effects of stress relaxation 9,27 (see Table I). In addition to S ad , both electrochemical and chemical water productions contribute to source term for ionomer water content. In other words, the produced water is assumed to be in absorbed (ionomer) phase. This is in agreement with the assumed structure for the CL in this work and has also been used by others. 42 Membrane gas crossover.-The model accounts for crossover of all gaseous species through the membrane. 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:  9 ( ,CL an mb where p i,CL an ca denotes partial pressure of gas species i (H 2 , O 2 , or N 2 ) in the anode/cathode CL, Ψ i denotes the corresponding membrane gas permeability, which is given in Table II for different species, and δ mb is the membrane thickness. Lastly, the nitrogen partial pressures can be approximated by their respective values in the flow channels with small errors. However, partial pressures of oxygen and hydrogen must be evaluated in the CLs due to their non-negligible fluxes through the thickness of the porous layers.
Reactant and water vapor transport in porous layers.-We only consider diffusive transport for the gas phase (reactant and water vapor) and neglect convective transport, since its effects are shown to be negligible under most typical operating conditions. 44,45 For reactant transport, we note that only concentrations at the active sites are of interest. These concentrations are determined based on the corresponding channel concentrations and the molar flux of each reactant, which is calculated from a known local current density.  where the Henry's law constant appears explicitly for H 2 , 46 but is implicitly included in  O2 for O 2 47 (see the subsection on Local transport resistance in the CL). The reactant transport resistance is calculated as the sum of convectional resistance through the flow channel and diffusional resistances through the different porous layers: The first term on the right hand side denotes the convective transport resistance at the GDL-CH interface, determined by the Sherwood number 48 : The second to fourth terms in Eq. 10 describe the diffusional gas transport resistance through different layers given by: where δ layer is the particular layer thickness and D i eff is the effective diffusivity of species i in that layer, which is calculated by accounting for molecular and Knudsen diffusion. Corrections are then applied to account for tortuous diffusion pathway through porous media. The molecular diffusivity values for different species are given in Table III and the effective diffusivity calculations are presented in Table IV. The last term in Eq. 10 is the local transport resistance in the CL. In this work, this term is assumed zero for all gas species except O 2 , whose local transport resistance is nonnegligible as discussed in the subsection on Local transport resistance in the CL.
Regarding water vapor transport, we again focus critically on the vapor concentrations in the CLs, as their values impact ionomer water uptake. Accordingly, the vapor dynamics in the CLs are given by: These equations account for vapor diffusion, water exchange with the ionomer, and the phase change process, where  v an ca is the total resistance to vapor transport between the CL and the channel (calculated by Eq. 10).
To capture the vapor distribution throughout the rest of the MEA, the quasi-steady state assumption is utilized due to the fast diffusive transport. More specifically, we initially disregard the phase change source term and obtain a linear distribution of vapor concentration through the thickness of all layers, namely, the CL, MPL, and GDL. This can be done by utilizing the diffusional resistances in each layer as well as the known boundary conditions at the CL and channel. In particular, we can formulate this as the following system of linear equations to be solved for vapor concentrations at the interfaces of various layers:   ( ) The above system is obtained by applying flux continuity conditions across different layers, i.e. the vapor flux through the CL should be the same as that through the MPL and GDL. The two known quantities in the above system are the boundary conditions, i.e. the vapor concentration at the flow channel and the CL (known from initial guess). The solution of the system will determine the vapor concentrations at the interface of the layers. This system is solved for both the anode and cathode sides. Overall, this allows us to obtain an approximate distribution of vapor concentration through the MEA thickness. Once this distribution is available, any supersaturated vapor results in liquid accumulation through condensation and the vapor fluxes are updated accordingly.
Regarding the source terms in Eqs. 13 and 14, S ad is the absorption/desorption source term described in Eq. 6 and S pc is the phase change source term given by: where k evp/cnd denotes the phase change rate, c sat (T) is the temperature dependent saturated vapor concentration, and c v is the vapor concentration. As the pore size becomes smaller, a larger interfacial area is available for phase change. 52 To model this effect, the phase change rate in the CLs is chosen to be the highest (50000 1/s), followed by that in the MPLs (1000 1/s) and GDLs (500 1/s). Moreover, high rates are used to ensure the facile phase change kinetics as suggested in the literature. 53 The dependency on the liquid saturation is used to model the fact that the interfacial area available for phase change scales almost linearly with liquid saturation. 53 Finally, 0.05 is used to capture the effects of a hydrophilic pore network that allows for phase change to occur at very low liquid saturation levels. 54 The saturated vapor concentration is defined as: where p sat is the temperature dependent saturation pressure 55 : where the solution variable is the liquid pressure (p l ). Assuming constant gas phase pressure in the through-the-membrane direction, the variations in liquid pressure uniquely determine the variations in the capillary pressure (p c ), which is defined as the difference between the liquid and gas phase pressures, i.e. p c = p l − p g and ∇p c = ∇p l . With known p c , the liquid saturation (s) can be determined. For this purpose, the Leverett J-function is used as the closure equation 56 to relate capillary pressure and liquid saturation: It is worth mentioning that the Leverett J-function was originally developed for homogeneous sand and its applicability to fuel cell porous layers is questionable. In light of such observations, more detailed models for the liquid water transport have been developed in recent years, the most notable of which is the mixed wettability model 52,57-59 that has become popular. 10 Nevertheless, recent literature also points to the fact that a representative elementary volume (REV) cannot be clearly defined in the through-plane direction for fuel cell porous layers. 60 Therefore, even the most advanced macrohomogeneous models do not have sufficient fidelity to capture the details of liquid water distribution in these layers. Accordingly, the Leverett J-function approach is adopted here for simplicity. Indeed, this model can be parameterized to capture the aggregate behavior of liquid accumulation in the porous layers with acceptable accuracy for real-time modeling purposes.
Heat transport.-The temperature distribution across different cell layers is captured with the heat equation: where k T eff is the effective thermal conductivity of the layers given in Appendix Table A·I. The membrane thermal conductivity depends on its hydration state and is given by 61 : The source terms for energy conservation are given in Table V. In Table V, H pc is the enthalpy of phase change 13 : Furthermore, H ad denotes the heat associated with ionomer water sorption/desorption 62 : Lastly, H HOR/ORR is the heat of reaction given by 58 : where χ is the cathode CL effective utilization factor (see the subsection on Proton and gas transport trade-off in the cathode CL).
Local transport resistance in the CL.-A major consideration for modeling transport phenomena in the CL is the local gas transport resistance. This resistance determines the efficiency with which the available reactant gas in the CL reaches the active catalyst site and has received significant attention in recent literature, 47,63,64 mainly due to efforts aimed at reducing the Pt loading. A more in-depth review of such efforts can be found elsewhere. 9 We only highlight the need for the model to capture such local transport resistance to be able to provide accurate estimates of variations in the performance due to changes in the CL structure. This is especially important for the cathode CL, as the local O 2 transport resistance is considerable. On the other hand, this resistance may be ignored on the anode side, since H 2 transport is usually efficient.
In light of the above discussion, we disregard the local H 2 transport resistance in the anode CL and only account for local O 2 transport resistance in the cathode CL. This is achieved through a homogeneous CL model originally developed by Hao et al., 47 which assumes the Pt particles to be deposited on carbon support and covered by an ionomer thin film. Condensation in the pores can result in a thin film of liquid water on top of the ionomer. This structural view is used to calculate the volume fraction of different phases in the CL, including Pt, carbon, ionomer, and pore space 47 : The ionomer and liquid water film thicknesses are given by: Following the above structural assumptions, the local O 2 transport resistance is assumed to be comprised of several components 47 : The right hand side terms in the above equation describe the resistance due to dissolution in water, diffusion through the water film, dissolution in ionomer, diffusion through the ionomer thin film, and adsorption onto the active Pt site, respectively. An important assumption of the model is that the resistances due to dissolution and adsorption are proportional to the diffusional resistances. 47 This assumption is made due to the lack of measured data for these interfacial resistances. In particular, three fitting parameters k 1 , k 2 , and k 3 are introduced 47 : Here we directly use the original formulation with the fitted parameter values from Ref. 47. However, these interfacial resistances may be lumped into a single effective resistance to simplify the fitting procedure. The reader is referred to the original publication 47 for the fitting procedure and further details about the homogeneous CL model.
Proton and gas transport trade-off in the cathode CL.-While the local O 2 transport resistance in the CL is captured with a homogeneous model, efficiently modeling the trade-off between proton and O 2 transport in the cathode CL remains a challenge. A distributed model of the CL that solves for the variations of O 2 concentration and ionomer conductivity through the thickness 9 will automatically explore such trade-offs. However, such an approach will be computationally intensive, as it requires keeping track of a relatively large number of states. Therefore, the challenge is in obtaining an approximate solution to this problem with the minimal number of state variables. Specifically, we seek a solution that requires only an average ionomer water content (according to Eq. 3) and liquid saturation in the cathode CL.
To obtain an efficient solution, we define a parameter for cathode CL utilization (χ). This parameter is used to explore the proton and gas transport trade-offs in the cathode CL. Particularly, χ represents the fraction of the cathode CL that is accessible to both O 2 and protons without significant mass or proton transport limitations. To calculate χ, one has to minimize the sum of ohmic and ORR activation losses in the CL, as shown schematically in Fig. 2. This is done by first estimating the effective proton transport resistance in the cathode CL using the method proposed by Neyerlin et al. 65 Shortly, the model calculates this effective protonic resistance as a function of the total resistance: Figure 2. Schematic for calculation of cathode CL utilization factor: The effective reaction front, shown with the vertical black line, is defined to be at the location of minimum voltage loss due to a combination of mass and proton transport limitations. It is assumed that no ORR takes place to the left of the effective reaction front due to severe mass transport limitations.  where R CL ca and R CL ca,eff denote the total and effective cathode CL protonic resistance, respectively, n Brugg is the Bruggeman exponent that is used to model the connectivity of the ionomer network, 33 and ζ is a polynomial fit to the solution obtained by Neyerlin et al., 65 given by: 3), which is expected to impact the estimated resistance to some extent. This difference is neglected here for simplicity.
If O 2 transport is efficient and imposes no limitations, i.e. χ = 1, then the above calculation yields the same estimate for the effective protonic resistance in the cathode CL as that obtained by Neyerlin et al. 65 In case of O 2 transport limitations, i.e. χ < 1, it is assumed that a fraction 1 − χ of the CL close to the membrane is not accessible to O 2 . Therefore, no reaction takes place in this region of the CL and protons coming through the membrane have to be transferred through this portion of the CL without participating in the ORR, i.e.
(1 − χ)δ CL is treated as an extension of the membrane. This results in an additional protonic resistance in the CL, which is taken into account through the second term on the right hand side of Eq. 36. Since the portion of CL closest to membrane does not participate in the ORR, the effective resistance to O 2 transport through the CL thickness has to be modified accordingly. In particular, we modify Eq. 12 for O 2 transport resistance in the CL as follows: As mentioned earlier, χ can be determined by minimizing the sum of activation and CL proton transport losses. This minimization can be carried out analytically as follows: where η ORR is the ORR activation overpotential (see Fig. 2). The solution is obtained through symbolic calculation, which yields a rather involved expression that is not presented here due to space considerations. The solution has to satisfy χ ∈ [0, 1]. In case the computed solution is out of this range for some conditions, it is projected onto this domain. For instance, χ > 1 yields the intuitive result of χ = 1. This completes the reduced model formulation for the cathode CL.
Terminal voltage and reaction kinetics.-The terminal cell voltage is given by: The expressions for different terms in the above equation are given in Table VI and the corresponding descriptions are provided in the Appendix Table A·I. For any given cell voltage, this equation is solved iteratively to find the corresponding current density at each location along the flow channel. The ORR kinetics model used in this work is based on the modified Tafel expression proposed by Subramanian et al. 66 that accounts for coverage of Pt sites with oxide species. The oxide coverage profile with respect to electrode potential as well as its dynamic growth are of critical importance in accurately modeling the transient performance of the cell. Accordingly, we use the following model for oxide growth dynamics, which assumes that water oxidizes Pt sites 12,67 : In the above equation k PtO is the oxide formation rate constant and η PtO is given by: where E rev,PtO is the reversible potential for PtO formation. Note that modeling the oxide film growth is especially important in capturing the voltage/current transients at higher potentials. Finally, the ORR kinetics are assumed to be RH dependent based on the experimental observations in the literature. 68,69 Although the significance of this effect has been questioned, 70,71 the approach adopted herein is still beneficial in capturing the impact of lowered RH on accessibility of Pt particles deposited inside the inner catalyst pores of high surface area carbon support. 72,73 This RH dependency is captured through the factor f RH that scales the cathode exchange current density. This factor depends on the ionomer water content as follows: where σ f is a sigmoid function defined as: In the above equations, f RH,a , f RH,b , and f RH,c are fitting parameters and λ sat (T) denotes the ionomer water content under saturated condition at any given temperature (see Table I). A similar approach has been suggested in recent PEM fuel cell modeling literature, wherein the ECSA is varied as a function of RH. 16,74 Flow channel model.-The flow channel model is a set of ODEs that determine the flow rates of reactant, nitrogen, and water in each channel. In particular, each ODE has the form of: where  M j denotes the molar flow rate of species j in the channel and the right hand side term (f j ) for each species is given in Table VII. In the table, the subscripts CH and LN correspond to the channel and land locations, respectively. In short, the channel model accounts for reactant consumption, gas crossover through the membrane (see the subsection on Membrane gas crossover, and water exchange with the MEA model.
To complete the channel model, the temperatures in both channels are assumed to be the same as the coolant temperature. Moreover, the coolant temperature is assumed to vary linearly from inlet to outlet based on the local current density. This is based on the fact that most of the heat generated in the cell is due to the reversible and irreversible heat of the ORR, whose dependence on current density can be approximated to be linear. The coolant temperature in segment i in the channel is given by: where T cool,in and ΔT cool denote the inlet coolant temperature and the change in the coolant temperature along the channel, both of which are controlled values and known, and N seg is the number of segments along the flow channel. In addition to the above assumption for the temperature variations, a linear pressure distribution is assumed from inlet to outlet in each channel, based on the knowledge of inlet and outlet pressures. This assumption is valid when slugs of liquid water do not exist in the channels, which is the case under most robust operating conditions used for automotive stacks.
Boundary conditions.-The heat equation is solved with convective boundary conditions at the coolant channel interfaces on both the anode and cathode sides. The convection coefficient (h conv ) is identified during the model parameterization step. For liquid water transport (Eq. 19), a zero flux boundary condition is used at the membrane interface with the CL, while the channel boundary condition is given by: Numerical implementation.-The PDEs for heat and liquid water transport (Eqs. 19 and 21) are spatially discretized using the control volume (CV) method. The through-plane direction is discretized using 20 CVs for each GDL, 4 CVs for each MPL, 4 CVs for each CL, and 1 CV for the membrane. The Backward Euler method is used for time stepping in both cases. For the liquid water transport problem, Jacobian linearization is utilized for the nonlinear terms. 75 The resulting linear system in each case (heat and liquid water transport problems) is tri-diagonal and is solved using the Thomas algorithm. 76 It should be pointed out that the relative permeability is zero when no liquid exists within a control volume. This eliminates the diffusive term and causes stability issues. In recent PEM fuel cell modeling literature, this issue has been overcome by adding a small artificial value to the permeability to ensure it never vanishes. 11,77 While this approach is certainly acceptable given that the added values are much smaller than the uncertainty in the actual permeability estimates, it is not a consistent stabilization method, since it increases the permeability artificially. In this work, we utilize permeability upstreaming as a consistent way of dealing with such instability issues. 78 Upstreaming essentially amounts to evaluating the relative permeability using the liquid saturation of the control volume in the upstream, which ensures that the permeability at the advancing front is never zero.
As for the water balance problem across the CCM, the set of ODEs given by Eqs. 2-4 and 13-14 are discretized in time using the Backward Euler method. This was deemed necessary as the equation system is stiff and the numerical damping offered by backward Euler helps in this regard. 79 The discretized system is solved along with the algebraic equation for the cell voltage (Eq. 40) using the Newton method with analytical Jacobian. Finally, forward Euler is used for the channel model (Eq. 45 and Table VII). Regarding the spatial discretization along the flow channels, either 12 or 27 nodes are used for the presented results depending on the channel length, as we have found that a channel segment that is larger than 1 cm in length can result in numerical issues in our implementation. Finally, a fixed time step size of 100 milliseconds is used for the model.
A version of the model has been implemented in C. The code is then used as an S-Function block in Simulink in order to simplify the user interface. The P2D 2 Simulink model with co-flow channel configuration and 12 nodes along the flow channel runs about 50 times faster than real time on a laptop computer with a 2.4 GHz processor. The worst case computation for a given time step typically takes less than 10 milliseconds. The counter-flow configuration is typically between 2 to 3 times slower than the co-flow case due to the iterative scheme employed (see below). All of these computation times are cut in half when the land region is not considered, i.e. a P2D model. This is the case for the model validation presented in Section 3. In addition, early-stage implementation with a MicroAutoBox has confirmed real-time simulation capability on dedicated hardware. Nonetheless, further work is needed to reduce the memory requirements for on-board implementation of the model.
Potentio-dynamic vs Galvano-dynamic simulations.-The model is formulated assuming potentio-dynamic mode of operation, i.e. cell voltage is the input to the model and current density is one of its predicted outputs. This approach allows us to obtain the current distribution along the flow channels without the need for costly iterations on the entire model. In particular, note that the electrodes must remain equipotential along the flow channels due to the high conductivity of the bipolar plates. Therefore, while the local current density varies along the flow channel depending on the conditions, the voltage remains constant. When voltage is considered as the  On the other hand, when current is the input to the model, the local current densities must be iteratively varied to ensure that the cell voltage remains constant along the flow channels. Therefore, modeling the galvano-dynamic mode of operation typically requires an iterative method, where the iterations are no longer restricted to a single segment (due to the global voltage constraint). Such an iterative scheme naturally results in a loss of computational efficiency.
To avoid this issue while simulating galvano-dynamic operation, we close the loop on the potentio-dynamic model using a simple proportional-integral (PI) controller to adaptively change the set point for cell voltage so that the predicted current density matches the galvano-dynamic input current density. In other words, the PI controller is used to find a suitable voltage input for the model. This approach, shown in Fig. 3, ensures zero steady-state error due to the integral action. Nevertheless, some transient error exists, which can be minimized by properly tuning the controller gains. In our simulations we have found that zero tracking error is achieved within 2 seconds of load change for all of the tested conditions. This time is in fact less than 1 second for most of the cases. Therefore, the transient loss of accuracy is minimal and the computational benefits clearly outweigh such losses.
Iterative scheme for counter-flow configuration.-The model also has the capability to simulate the counter-flow configuration. This is achieved by designing a new iterative scheme. In particular, the anode stream flows in the opposite direction of the cathode and coolant streams under the counter-flow configuration. Therefore, this case cannot be handled easily within the P2D modeling framework, as this framework relies heavily on the flow of information on both sides of the cell being in the same direction. This necessitates the need for an iterative scheme. We have designed the iterations around the total water flow rate in each segment of the anode channel: we start with an estimate of the vapor flux at the channel-GDL interface at each location along the anode channel. These flux values are used to calculate the total water vapor flow rate and the corresponding vapor concentration in each channel segment. The vapor concentration values in the anode channel can be used to solve the model similar to the co-flow case and calculate the vapor fluxes at the anode channel-GDL interfaces. The calculated values are compared with the initial guessed values and the process is repeated until convergence is achieved. This iterative scheme is shown in Fig. 4 and is numerically stable and converges (with a relative tolerance of 2%) in less than 10 iterations in most cases. With this number of iterations, the model can still achieve real-time performance even when the worst case computations are considered.

Model Validation
Experimental data.-The model is validated against data from two stacks with different material sets and designs. The detailed operating conditions and stack designs are proprietary. However, a high level description of each stack and the corresponding dataset is provided below.
The first dataset, denoted by "Stack A", contains polarization and high frequency resistance (HFR) data (measured at 1 kHz) obtained on a cell with a 0.15 mg· cm −2 cathode Pt loading and a small active area (40 cm 2 ) under differential operating conditions (high stoichiometric ratio). This dataset, which is collected under co-flow configuration, includes 4 different operating conditions intended to replicate various real-world scenarios and constitutes a rather wide operating range, covering temperatures from 40 to 80°C. Importantly, high current densities up to 3 A· cm −2 are achieved in each case, allowing for the model predictions under such high loads to be critically examined. Due to differential operation setting, no significant transients are observed in the data. Therefore, model validation is carried out using the steady state polarization and HFR values from Stack A.
The second dataset, denoted by "Stack B", contains only polarization data from an automotive size stack with a 0.4 mg· cm −2 cathode Pt loading and a large active area (about 300 cm 2 ). It should be noted that the material set (membrane, porous layers, catalyst loading, and CL morphology) used in Stack B is considerably different from that of Stack A. Additionally, different operating conditions are used, where the effects of humidity conditions have been extensively explored at 70°C. Moreover, variations in the stoichiometric ratios and low stoichiometric operation under integral conditions allow us to better investigate mass transport limitations compared to Stack A data. Each loading condition on the polarization curve is maintained for a sufficiently long time to establish steady state operation before the load is changed in a step manner.
Note that the plates used in both stacks have very small land/ channel ratio on both the anode and cathode sides. As such, the inplane distributions are disregarded in the model, i.e. the transport equations are only solved under the channel (P2D model). The role of transport under the land is explored with the P2D 2 model in the Simulation Case Studies section.
Model parameterization.-A set of 20 parameters are identified using experimental data for each stack. This particular choice was made based on a preliminary sensitivity analysis that showed the   Table VIII along with their values for each stack. The parameter descriptions can be found in the appendix. As seen in Table VIII, the parameters are constrained within specific ranges that are based on the uncertainty of the corresponding values reported in the literature. These constraints help avoid over-fitting the model to data and ensure the identified parameters remain physically meaningful. Other parameters of the model are provided in the Appendix Table A·I.
The parameters specified in Table VIII represent factors that control the kinetic, ohmic, and mass transport characteristics of the cell response. Therefore, they can be identified separately using data from different parts of the polarization curve. 80,81 However, in this work we identify all of the parameters together using the entire dataset. In particular, we seek parameter values that minimize the following cost function: where ⊙ denotes element-wise multiplication, e V is the vector containing n V voltage predictions errors in mV, e R is the vector of n R HFR prediction error (for Stack A data) in mΩ · cm 2 , r 1 and r 2 are weight parameters (both set to 100 in this work), and I(x) denotes the indicator function of x, given by: With this, the last two terms in Eq. 48 put an additional penalty on voltage errors exceeding 10 mV or HFR errors exceeding 10 mΩ · cm 2 . This cost is minimized for each stack separately, identifying a different parameter set for each stack. Note that the steady state data is used for parameterization. Since the model is transient in nature, it is allowed enough time to reach the steady state conditions before comparison is made with the experimental data. Moreover, since the parameters have largely different ranges, they are all scaled to a unit interval (between 0 and 1) using linear or logarithmic minmax scaling (see Table VIII): whereq i is the normalized parameter value and θ i,min/max denotes the respective parameter bound. This scaling was determined to be essential in successfully identifying all of the parameters together.
The identification is carried out using particle swarm optimization (PSO) 82 in MATLAB, which is a global optimization routine that efficiently explores the entire parameter space to find the optimal parameter set. A swarm size of 200 was used and the optimization was terminated after 200 iterations, at which point the cost was not decreasing any further. The algorithm was executed several times for each dataset and the results with the lowest cost are reported here. It should be noted that due to the stochastic nature of the algorithm, each run can yield a new parameter set.
Results for stack A.-The polarization curves for the 4 conditions comparing the experimental data and the model predictions are shown in Fig. 5. The corresponding results for the HFR are provided in Fig. 6. Each polarization curve consists of 15 points, starting from OCV and ending at 3 A· cm −2 . This results in a total of 60 polarization data points across all 4 operating conditions. As for HFR data, only 7 measurements were made, shown as diamonds in Fig. 6. The experimental trend lines shown in the figure are estimations based on prior measurements using the same test stand.
Overall, we observe that the model captures the trends in the data very well. The maximum error in voltage across all 60 data points is 11.9 mV, while the maximum HFR error is 10 mΩ · cm 2 . Such a good prediction capability is partly due to the extensive model parameterization effort. Nonetheless, the fact that these fits are obtained with reasonable parameter values (see Table VIII) illustrates the model's capability in capturing the most relevant physical phenomena critical to performance predictions. This indicates that the model is indeed capable of describing the physical phenomena under the operating window for which it is parameterized.
A few comments about the parameterization process are in order. First, we found that the polarization data alone may be fitted with a variety of parameter combinations. In other words, the model is not structurally identifiable 83 using only voltage measurements at these 4 operating conditions. The fact that a polarization curve may be reproduced with different parameter combinations has also been pointed out by others. 84,85 Nonetheless, the operating conditions can be optimized to ensure maximum information is contained in the polarization data. 80 In fact, even in this dataset, having the polarization data at 40 and 80°C helps reduce the uncertainty in the identified parameters considerably. Second, as pointed out recently by Vetter et al., 34,86 most of the consequential parameters in these types of PEM fuel cell models are related to the membrane water uptake and transport. Therefore, adding the HFR data to the identification process can considerably increase the confidence in parameter estimates. These issues are extensively studied in separate publications 81,87 and are not pursued here any further.
Results for stack B.-The validation results using data obtained from Stack B with co-flow configuration are shown in Fig. 7. Overall, 5 different humidity conditions were tested, ranging from 30% to 100% RH, where the RH value indicates humidity conditions in both the anode and cathode channels. The data at three RH conditions were used in the optimization routine to parameterize the model. The same parameter values were then used to validate the model predictions against experimental values for the other two humidity conditions. As seen in Fig. 7, the model predictions match the experimental results closely. The largest error in voltage prediction is 12.8 mV. This indicates the model's capability in accurately predicting the stack performance when parameterized properly. Nevertheless, we note that during model parameterization, we were able to find several parameter combinations that yield reasonable model predictions as compared to available voltage measurements. This further signifies the structural identifiability issues that are common for physics-based fuel cell models and are studied in separate publications. 81,87 The transient behavior for the dry conditions (RH = 30%) is shown in Fig. 8. The figure depicts voltage dynamics from experimental data with co-flow configuration along with the corresponding model predictions. The model captures most of the features seen in the transient measurements. In particular, the nonmonotonic voltage response (seen as overshoot and undershoot) to step changes in the current density is mostly captured. This response is due to membrane hydration/dehydration at higher current densities, while Pt oxide coverage dynamics is the main contributing factor at lower current densities. Note that the model predictions for the immediate response to current density variations are also affected by the PI controller performance. Therefore, some discrepancies between the model and measured transients immediately following a step change in current density are to be expected. Since the transient data were not used during parameter identification, the good agreement between the transient model predictions and the corresponding measurements serves as a further validation of our parameterization approach.
In Fig. 8, two sections of the voltage transient response are magnified and provided in the insets for further discussion. The first part corresponds to the highest load used in the experiments (1.4 A· cm −2 ). During this time, some voltage swings are observed, which correspond to variations in the inlet coolant temperature. A decrease in the inlet coolant temperature is observed to result in an increase in the cell voltage, which is attributable to better membrane hydration. The voltage sensitivity to these temperature variations is accurately captured by the model. The second inset shows the voltage dynamics at lower loads above 0.8 V. For this case we observe a voltage overshoot followed by a relatively slow decay in the measured cell voltage after the current density is lowered. This behavior is mainly due to the formation of oxide species on Pt sites and is partially captured by the model (see the Pt oxide coverage plot in Fig. 8). In particular, the measurements show signs of a second-order system with two separate time constants; a fast decay on the order of tens of seconds, followed by slower response that settles on the order of hundreds of seconds. This behavior may be attributable to formation of different oxide species. 88 Since the present model only incorporates PtO coverage, the slow dynamics are not captured accurately. Nonetheless, such behavior rarely has relevance to performance optimization for it to prompt further model development.
In addition to voltage dynamics, Fig. 8 also shows the membrane hydration, average membrane temperature, and current distribution along the flow channel, as well as the oxide coverage dynamics. Membrane hydration, temperature, and current density are all seen to increase along the flow channels. However, membrane hydration shows an interesting behavior, where it dramatically increases near the outlet region. This is especially the case for drier conditions with the co-flow configuration. Under such conditions, the low inlet humidity results in a dry membrane. But the water produced near the  inlet of the cell humidifies the channel gas streams, which in turn helps improve hydration near the outlet region.
To finalize our discussion on the modeling results for the experimental conditions, Fig. 9 provides the along-the-channel distributions of current density, membrane hydration, normalized water crossover, and channel and CL RH at a current density of 1.4 A· cm −2 (corresponding to Time = 1700 s in Fig. 8). The plots also provide a comparison between distributions for the co-flow and counter-flow configurations. A few observations are in order. First, we note the close correspondence between the current density and  membrane hydration distributions along the flow channels. Particularly, the membrane hydration is seen to increase monotonically from inlet to outlet for the co-flow case, whereas a peak in membrane hydration is observed near the middle of the cell for the counter-flow configuration. As for the current densities, the distribution for the counter-flow case aligns with its corresponding membrane hydration. For the co-flow case, the current density is seen to increase from the inlet toward the outlet with improved membrane hydration. However, near the outlet, the current density drops slightly due to mass transport limitations. In fact, this current distribution is typical for co-flow operations and the location of the peak current density is closely related to the inlet humidity conditions, with higher humidities resulting in a peak current density closer to the inlet. The figure also shows the normalized membrane water crossover defined as the rate of water crossover to the anode side normalized by the local water generation: As is typically the convention, a positive value denotes flow from anode to cathode. This plot clearly shows the water recirculation mechanism for the counter-flow configuration, wherein near the cathode inlet the humidified anode outlet stream results in a net water crossover toward the cathode. However, as we move along the cathode channel and the cathode stream becomes humidified, the net water crossover changes toward the anode, humidifying the incoming dry anode stream. This mechanism, which can also be seen from the channel RH distributions, is the main feature of counterflow operation reported in both experimental and modeling literature. 18,89,90 However, the presented model is the first model to capture this behavior in a computationally efficient manner. The same plot also shows that the net water crossover is consistently toward the anode side for the co-flow case, which is in agreement with experimental observations for state-of-the-art and thin membranes. 91 Finally, the CL RH values are also shown in Fig. 9. This plot highlights the fact that even when the channel streams are saturated, the CL can experience sub-saturated conditions due to the heat released through the electrochemical reactions. This has also been observed experimentally in the literature. 92 While measurements for such distributions were not available in our experimental set-up, their qualitative agreement with experimental results in the literature showcase the present model's capabilities in capturing these trends. Overall, the very good quantitative agreement with available measurements, the aforementioned qualitative agreement with literature data, the robustness of the numerical algorithms in solving the model for real operating conditions without convergence issues, and the fact that these results are obtained more than an order of magnitude faster than real time all point to the utility of the present model and the underlying numerical approach to act as a software sensor for PEM fuel cell systems.

Simulation Case Studies
With the model validated with various experimental data in the previous section, we now turn our attention to some of its prediction capabilities using some load cycling scenarios. Note that while a P2D version of the model was used for validation, the P2D 2 model is used in this section to demonstrate the ability of the bi-domain modeling approach in capturing the in-plane distributions. Here, the current density is cycled between low and high loads under two distinct conditions. Simulations for each condition are conducted with both co-flow and counter-flow configurations. The operating conditions are shown in Table IX. When the anode and cathode conditions are not separately identified (e.g. RH), they are taken to be the same. Each loading condition is maintained for 300 s before a step change in the load. This time allows for the cell dynamics to settle to a quasi steady state. Moreover, note that the step change in the load is chosen to be large enough to excite the system dynamics. Similarly, the two sets of operating variables (cold wet and hot dry) are chosen to show the model's prediction capabilities under extreme conditions. A complete list of parameter values used for these simulations is provided in the appendix.
The average dynamic responses of the cell under the various simulation scenarios are shown in Fig. 10. As seen in the figure, the step increase in the load results in a very small undershoot in voltage, which is then followed by a slow decay. The initial undershoot can be attributed to a combination of EOD that dries out some parts of the membrane as well as the existence of the PtO layer that has accumulated at higher voltage. The slow decay in voltage that follows is seen to be a direct consequence of membrane dry-out. This can in turn be attributed to the significant heat that is generated at this higher load, as well as the high EOD. The figure also shows that the flow configuration does not have a significant impact under cold and wet conditions, whereas counter-flow configuration considerably improves performance under hot and dry conditions. This is to be expected, since counter-flow operation yields better membrane hydration conditions as can also be seen in the figure. Figure 10 further shows the average cathode CL utilization (χ) for the different scenarios. It is observed that the cold and wet conditions result in lowered CL utilization. This is due to mass transport losses incurred by the accumulated liquid water under these conditions. The CL is also not fully utilized under the high load for the hot and dry conditions with co-flow configuration. Although the average utilization is very close to unity, this amounts to a very interesting observation, namely, dry ionomer in the cathode CL can also cause local mass transport losses and reduce the catalyst available for ORR. Of course, this behavior can already be captured by comprehensive fully coupled models. 9 Nevertheless, computationally efficient solutions with this level of fidelity are not commonly pursued in the literature. Figure 10 also provides the distinct current density profiles under the channel and land regions for two of the simulated cases. As seen in the figure, the current density is lower under the land region for the cold and wet condition, which is due to the mass transport limitations. However, for the hot and dry conditions the current production under the land is on par with that under the channel, since better hydration under the land compensates for increased mass transport resistance. Lastly, the figure shows liquid saturation in the cathode CL and GDL under the land and channel regions for the cold and wet conditions. Clearly, the land region is more susceptible to flooding due to lower temperatures and higher mass transport resistance. This is in qualitative agreement with more complex modeling efforts. 9,54 More interesting, however, is the dynamics of liquid accumulation in the CL and GDL. The CL shows a tendency to flood rapidly, within 20 seconds of the load change, whereas liquid accumulation in the GDL is seen to be much slower, not reaching steady state under the land even after 300 seconds. These disparate time scales are due to differences in porosities of the layers, as well as the environmental conditions experienced by each layer.
To further demonstrate the model's prediction capabilities, example distributions obtained from the model for two of the simulated scenarios are shown in Fig. 11. The model is seen to capture the main features of the in-plane distributions, namely, higher mass transport resistance under the lands, which leads to increased flooding and improved membrane hydration, and the higher temperature under the channels. More specifically, the distributions for cold and wet conditions with counter-flow configuration show that liquid accumulates mostly under the land region. In this region, each GDL is more prone to flooding near the outlet of its respective gas stream. However, the cathode CL is mostly flooded near the cathode inlet, which is due to the significant water crossover from the anode side in this region. The good hydration on the anode side also results in improved membrane hydration near the cathode inlet. The figure also shows the along-the-channel current distributions for both the channel and land regions, where it is observed that current generation under the land is suppressed due to mass transport limitations. Another interesting observation relates to cathode CL utilization under the land and the fact that it is reduced near the cathode inlet due to cathode CL flooding in this region. Moreover, we observe that this reduced utilization has also resulted in higher temperatures in this region due to higher volumetric heat generation.
Lastly, Fig. 11b shows the same distributions for the hot and dry conditions under co-flow configuration. We observe that no liquid has accumulated anywhere in the cell. As for the current distributions under the channel and land regions, they go through a peak as we move along the flow channel. Near the channel inlet, the land region produces more current due to better hydration. Once the peak current is achieved, the current under the land region drops more sharply due to more severe mass transport limitations. These transport limitations also result in reduced cathode CL utilization under the land near the channel outlet.
These results demonstrate the prediction capabilities of the model presented here. As mentioned earlier, it is the balance between the fidelity and the low computational cost of the present modeling framework rather than the specific simulation results that is of significance. The fact that the model is able to reproduce the behavior that is captured by much more computationally intensive models verifies the approach adopted herein. Finally, the high fidelity and efficiency of the model offers significant advantage for real-time estimation as well as the tasks that require repetitive model evaluations such as optimization of operating conditions.

Summary and Conclusions
A high fidelity and computationally efficient model of PEM fuel cells is presented. The model takes advantage of the disparate time and length scales to spatio-temporally decouple the problem and achieve a computational time that is more than an order of magnitude faster than real time. Particularly, the model efficiently captures the in-plane and through-plane distributions utilizing a pseudo-2D bidomain modeling approach. Additionally, a novel reduced-order model for water balance across the MEA is presented and a simple analytical method is used to model the proton and mass transport trade-offs in the cathode CL and capture its effective utilization. These novelties further contribute to the efficiency of the presented model. An optimization-based approach is used to parameterize the model for two different fuel cell stacks and some issues regarding the model parameterization process are discussed. The model is validated against a variety of operating conditions for the two stacks and shown to capture the experimental performance measurements. Finally, simulation case studies are presented to showcase the model's prediction capabilities under different operating scenarios. Based on these results and initial implementation on dedicated rapid-prototyping hardware, it is concluded that the fidelity and computational

Acknowledgments
Financial and technical support by Ford Motor Company is gratefully acknowledged.

Appendix
A complete list of model variables and parameters is provided in Table A·I below for accessibility. We note that the parameter values shown here correspond to the parameters used in the Simulation Case Studies section. The anode and cathode parameters are assumed to be the same unless noted otherwise.