A computationally efficient and high-fidelity 1D steady-state performance model for PEM fuel cells

The performance of a proton exchange membrane (PEM) fuel cell is determined by many factors, including operating conditions, component specifications, and system design, making it challenging to predict its performance over a wide range of operating conditions. Existing fuel cell models can be complex and computationally demanding or may be over-simplified by neglecting many transport phenomena. Therefore, a high-fidelity and computationally efficient model is urgently needed for the model-based control of fuel cells. In this study, semi-implicit multi-physics numerical models have been established, taking the mass, momentum, reactants, liquid water, membrane water, electrons, ions, and energy in all fuel cell components into account. The developed 1D model is of high fidelity by incorporating the two-phase flow, non-isothermal effect, and convection, and is still computationally efficient. These models are validated against data from an auto manufacturer with good agreements, and the computing efficiency is evaluated on a modest laptop computer. The modeling results suggest that the two-phase flow model exhibits better prediction accuracy than the single-phase flow model when reactants are fully humidified, while under low humidity conditions, the two models present equivalent performance as liquid water does not exist in the fuel cell components. The results also suggest that the maximum convective/diffusive ratio of H2, O2, and vapor mass fluxes can be 12%, 5.3%, and 35%, respectively, which are ignored in most diffusion-dominant models. The developed models are computationally efficient, requiring only 0.56 s and 0.26 s to simulate a steady-state operation of fuel cells for the two- and single-phase flow models, respectively. This implies that the developed models are suitable for the control of PEM fuel cells.


Introduction
The proton exchange membrane (PEM) fuel cell has become an emerging clean energy technology that utilizes renewable fuels and exhausts zero onboard gas emissions for automotive applications [1]. The fuel cell performance can be affected by many factors, including operating conditions (e.g. relative humidity, temperature, pressure, flow rates, current density, and voltage) [2,3], component development (e.g. materials, compositions, structure, and dimensions) [4][5][6][7][8][9][10][11], and system design (e.g. reactant supply, water and heat management, and electric load) [12][13][14][15][16]. Therefore, the fuel cell operation is complicated with many degrees of freedom, and the reliable and stable operation of fuel cells requires a proper control strategy for industrial applications. To maintain the high operating performance, model predictive control is a prevalent control strategy in comparison with traditional proportional integral derivative control and adaptive control [17], which utilizes a fuel cell model to predict the output variables. Therefore, a control-oriented fuel cell performance model with high prediction accuracy and low computing cost is urgently needed in industry to implement model predictive control strategies and to avoid harmful operating conditions. The most prevalent models in the model-based control of fuel cells are 0D analytical and empirical models [18][19][20][21][22]. For instance, Tirnovan and Giurgea [20] developed an empirical fuel cell model based on the first law of thermodynamics, taking the pressure drop in pipes and channels into account. An electrochemical model is employed to adjust the flow rates of air to achieve maximal system energy efficiency. This model is computationally efficient due to its semi-empirical features without the gas and water dynamics. Danzer et al [23] developed a system-level fuel cell model based on a flow resistance network to prevent oxygen starvation by regulating the pressure, current density, and stoichiometry ratio. However, the fuel cell model validation was not reported, and its accuracy was unclear. Matraji et al [24] developed a nonlinear dynamic fuel cell model based on mass conservation and ideal gas law. The model is used to control the pressure of gaseous reactants and prevent component damage, by minimizing the pressure difference between anode and cathode with a multi-input multi-objective (MIMO) controller. Hong et al [25] established a MIMO nonlinear model based on simple ideal gas law and mass balance to optimize the reactant flow rates for PEM fuel cells with recirculation of hydrogen. As can be seen, the models widely used in the model-based control of PEM fuel cells are often analytical and empirical without considering component dimensions and transport behaviors for fast computing speed. However, these models, often based on conservation laws and state equations, were over-simplified, and factors affecting the PEM fuel cell performance, including convection, reactant supply, and water (liquid, vapor, and dissolved) dynamics, were neglected.
Many efforts have been devoted to accurately predicting cell performance by taking various transport phenomena and actual cell dimensions into account [26][27][28]. Many 3D models have been developed to understand the various factors affecting the overall performance of PEM fuel cells [2,26]. For instance, Zhao and Li (2019) developed a 3D single-channel model to investigate the oxygen transport mechanisms in fuel cell components, and the simulation results suggest that the oxygen transfer from a straight channel to reaction sites is dominated by diffusion while convection accounts for 5%-16% of the total oxygen transport [26]. Various physical and electrochemical phenomena have been modeled, including activation and ohmic losses, multispecies transport, water transfer, and thermal effects, which are difficult and expensive to be experimentally studied. However, the 3D model is usually computationally inefficient and takes from a few hours to months to simulate a steady-state scenario [29][30][31], making it unsuitable for control applications. Therefore, many reduced dimensional models, typically pseudo-2D and 1D models, have been developed to capture the nonlinear transport-related fuel cell behaviors (in comparison with 0D analytical and empirical models) and to enhance the computational speed (in comparison with 3D models) [12,[32][33][34][35][36]. Pseudo-2D and 1D models can take electrochemical kinetics, ohmic polarization, and transport of chemical reactants and products into account, although the computational speed is often not reported in these studies. Recently, Goshtasbi et al [32,37] built a semi-empirical, unsteady pseudo-2D model by treating catalyst layers (CLs) as dimensionless interfaces, omitting the transport phenomena and phase changes in CLs, and ignoring gas dynamics in all porous components. The simplifications accelerated the computation speed, which was about two times faster than the actual physical processes. The authors pointed out that the model is only valid for large time steps (e.g. ∼100 ms), where gas dynamics are omitted. Jiang et al [1] performed a parametric analysis of PEM fuel cells based on a 1D fuel cell model with assumptions of instant phase change between liquid water and vapor and uniform pressure in different cell components. A switching function is defined to predict the liquid water and vapor in the computational domain. This simplification may significantly accelerate the computational speed; however, its accuracy in predicting the water phase transfer rates and vapor states requires further validation.
It is seen that most existing 1D models ignore some important phenomena such as two-phase flow, non-isothermal effect, and convection [17], while these simplifications may cause inaccurate prediction of fuel cell performance under certain conditions. For instance, single-phase flow models may cause errors when the fuel cells are operated at high relative humidity and high current density conditions [38][39][40]. Yu et al [41] developed an single-phase model to control the preheating process of the cold-start of the PEM fuel cells, in which phase changes are ignored. The model is demonstrated suitable for cold start in dry conditions, while its accuracy in wet conditions remains unclear. Isothermal models may have errors in predicting the phase changes within the porous media as the local temperature variation directly impact the local saturated vapor partial pressure and the phase change rates [42,43]. Convection has a significant impact on the transport phenomena within the fuel cell porous components and may affect the prediction of concentration polarization and water vapor at channel outlets, which is omitted by most 1D models [1,39,40,43]. Gong et al [44] optimized a 150 kW fuel cell stack and a centrifugal air compressor based on a 1D non-isothermal fuel cell model through an adaptive optimization matching method. The phase change rate between liquid water and vapor is assumed to be gigantic although the phase change equations are not provided, and convection effect is ignored. The phase change equations are provided in Gong et al's [45] 1 + 1D model, in which the explicit scheme is implemented in C language. The diffusion and convection are ignored in fuel cell electrodes, which significantly simplified model and accelerated the computing speed. The authors mentioned that the computing speed is faster than real-time physical phenomena but is still limited by the tiny time step applied to avoid instability issues of explicit numerical schemes. Based on the literature survey, the phenomena of two-phase flow, non-isothermal effect, and convection have not been simultaneously and comprehensively incorporated in a 1D fuel cell model yet, and the computing speed of 1D models needs to be improved for control applications. It should be also noted that many data-driven models [46][47][48] have been developed for the fuel cell control and do not require any physical knowledge of the fuel cell stack. However, data-driven models require a high-volume and high-quality experimental dataset and can be fuel cell specific, which is out of the scope of the present study.
Therefore, the objective of the present study is to investigate the nonlinear behaviors of the PEM fuel cells under steady operation by incorporating the key transport phenomena, including but not limited to two-phase flow, non-isothermal, and convection, as well as electrochemical kinetics within all corresponding components in a 1D numerical model. To achieve a high computational speed, a semi-implicit method for pressure-linked (SIMPLE) algorithm is applied to solve the governing equations simultaneously at each iteration step and a hybrid numerical discretization scheme is applied to the space-dependent terms [49]. The developed model is validated against data from Toyota for Mirai vehicles, and the computing speed has been evaluated for different convergence criteria and mesh numbers.

Development of 1D numerical model
The 1D PEM fuel cell performance model proposed in this study considers the transport phenomena along the direction from anode flow channel to membrane (hereafter it is called the through-plane direction). The computation domain is a typical 11-layer PEM fuel cell stack, including a membrane, 2 CLs, 2 microporous layers (MPLs), 2 gas diffusion layers (GDLs), 2 straight flow channels, and 2 bipolar plates (BPs), as shown in figure 1. The key parameters of thickness, porosity, and permeability of major fuel cell components are presented in table 1.
The governing equations for the 1D steady-state fuel cell operation are summarized in table 2, with the source terms of each governing equation presented in table 3. The mass conservation is shown in equation (1) (see table 2), and the gas-phase mass consumption and generation due to phase change and electrochemical reaction are incorporated in the source term, S m (see table 3). The momentum equation is shown in equation (2), and the effect of pressure drop due to permeation in porous media governed by Darcy's law is considered in the source term, S u . The gas species (including hydrogen, oxygen, nitrogen, and water vapor) transported through the channels and porous electrodes are governed by equation (3) via convection and diffusion, and the consumption of H 2 and O 2 and generation of water vapor due to phase change are described in the source term, S i . The transport of liquid water is governed by the gas-liquid interaction and capillary pressure gradient via equation (4), and the liquid water generation or consumption determined by phase change is mathematically determined by the source term, S l . The transport of dissolved water in the ionomers of membrane and CLs is dominated by the electro-osmotic drag effect and diffusion (see equation (5)), and the dissolved water change is from either cathode half-cell reaction or phase change with water vapor as depicted in the source term, S d . The transport of protons and electrons is governed by equations (6) and (7), respectively, via Ohm's law of conduction. The energy conservation is determined by equation (8), which is transported via conduction coupled with heat sources (S T ) of Joule heating, reaction heat, and phase changes in the corresponding components. A more detailed explanation of the governing equations for the fuel cell models employed in this study can be found elsewhere [17].
The fuel cell output voltage is calculated as [2,50], where V cell is the output voltage in (V); η act is the activation overvoltage in (V) within anode and cathode CLs; and η ohm is the total ohmic overvoltage in (V); V rev , is the thermodynamically reversible voltage in (V), determined by its reactants' and products' conditions and properties [51][52][53] via Nernst equation:  where g represents the Gibbs function in (J kmol −1 ), s denotes entropy in (J kmol −1 K −1 ), and the subscript, ref, represents the properties evaluated at the standard temperature (25 • C) and pressure (1 atm). The rates of O 2 and H 2 consumption, H 2 O generation, and heat production are determined by the electrochemical reaction rate in the CLs. The reaction rate is governed by the CL structure, reactant Gas species Dissolved water Note: GDL, MPL and CL denote gas diffusion layer, microporous layer, and catalyst layer, respectively; S is source term in different governing equations; concentration (H 2 in anode and O 2 in cathode CLs), liquid water saturation, activation polarization, and exchange current density. The reaction rate is often expressed by the Butler-Volmer equation [54,55]. To reduce the computational load, the Butler-Volmer equation is often simplified under various assumptions when it is applied to numerical models [56]. The hydrogen oxidation reaction (HOR) kinetic is often linearized based on the assumption that the HOR is fast with a small surface overpotential, while the oxygen reduction reaction is sluggish and hence of high overpotential, which allows the neglect of the anodic reaction term of Butler-Volmer equation [56]. Therefore, the Butler-Volmer equation may be reduced to the following form where α is the transfer coefficient, and R is the universal gas constant in (J kmol −1 K −1 ). The major physical and electrochemical properties and corresponding empirical correlations employed in the model are given in table 4.

Numerical procedure
The numerical model is implemented via Matlab by discretizing governing equations in every cell component and solving the discretized governing equations simultaneously based on the SIMPLE algorithm [49]. The governing equations are discretized by a hybrid central and upwind differencing scheme (UDS). The central differencing scheme (CDS) is employed when the diffusion dominates in the mass transport for small Peclet number (Pe < 2), while the UDS is adopted when convection dominates for large Peclet number (Pe > 2). The CDS solution is accurate to the second order of the Taylor series truncation error but may exhibit severe oscillations with coarse meshes, and the UDS solution is accurate to the first order but accounts for the flow direction. Therefore, a hybrid differencing scheme is employed in this study. The SIMPLE algorithm follows a guess-and-correct procedure for the 1D steady-state two-phase flow model, as shown in figure 2. The pressure, velocities, and other variables are initialized with guessed values, and then momentum equations are solved to correct the pressure and velocities. The other governing equations including mass, gas species,  water-vapor-to-liquid, dissolved-water-to-vapor, and dissolved-water-to-liquid, respectively; and the subscripts of a, c, g, ele, ion, eff and act stand for anode, cathode, gas, electronic, ionic, effective, and activation, respectively.  Note: T is the temperature in (K); p is pressure in (Pa); ε is porosity; Φ is volume fraction; K0 is the intrinsic permeability in (m 2 ); D i is the bulk diffusivity of species i in (m 2 s −1 ); µ is the dynamic viscosity in (N s m −2 ); λ is dissolved water content in ionomer in (kmol(H2O)·kmol(SO3H) −1 ); ω is the volume fraction of ionomer in catalyst layer; κ is the electrical conductivity in (S m −1 ); ∆s is the entropy change of reaction in (J kmol −1 K −1 ); j is the volumetric reaction rates in (A m −3 ); and the subscripts of a, c, d, g, l, v, ele, ion, ref and eff stand for anode, cathode, dissolved, gas, liquid, vapor, electronic, ionic, reference, and effective, respectively. liquid water, dissolved water, electrons, ions, and energy are simultaneously solved, and finally, the Nernst and Butler-Volmer equations are computed. If the solution is not converged, the variables will be corrected based on the results obtained at the previous iteration step. If the solution is converged, the calculation ends.
To understand the two-phase flow in the PEM fuel cells, a single-phase flow model is also established for better comparison. In the single-phase flow model, liquid water saturation in all fuel cell components is set to zero, meaning that the phase change and transport of liquid water are not computed.

Boundary and initial conditions
The boundary conditions are set based on the operational conditions, including cell temperature, inlet gas pressure, temperature, and flow rates, gas species composition and concentration, and electric loads. The cell temperature of the BPs is assumed constant, which is allowed to vary from 328-358 K from case to case depending on the available experimental data, and the inlet gas reactant temperature is identical to the BPs, assuming that the reaction heat is promptly removed by an efficient cooling system. The inlet pressure for the anode and cathode is directly based on the data obtained from Toyota. The mass flow rate of the reactants at the inlets is determined by the stoichiometry ratio. The molar concentration of the fuels and oxidants at the anode and cathode channel inlets, respectively, are calculated based on ideal gas law. The liquid water in the channels is assumed to be effectively removed upon its formation. Figure 2. The numerical procedure of the 1D steady-state PEM fuel cell model based on SIMPLE algorithm (p is pressure; u is superficial velocity; Φ represents other variables including mass, gas species concentration, liquid water saturation, the dissolved water content in ionomer, voltage, and other physical properties; and * represents the previous iteration step).

Convergence criteria and mesh independence
To evaluate the steady-state performance of the PEM fuel cells, the convergence criterion for the steady-state calculation is based on the maximum relative change of variables within all computational domains between two consecutive iteration steps via the following equation, where θ * and θ represent the variables at the previous and current iteration steps, respectively (e.g. voltage, temperature, pressure, liquid water saturation, dissolved water content, and species concentration). A parametric study to determine the optimized convergence criteria is performed and the results are compared with the data provided by Toyota (a random case is selected as a representative: current density of 2.0 A cm −2 , cathode oxygen stoichiometric ratio of 1.5, cathode inlet air pressure of 2.2 bar, anode hydrogen stoichiometric ratio of 1.25, anode inlet relative humidity of 100%, anode inlet pressure of 2.8 bar, and operating temperature of 328 K). For better comparison, the calculation errors in different variables of voltage, temperature, pressure, dissolved water content, liquid water saturation, and concentrations of oxygen, hydrogen, water vapor, and total species in various computing domains with fine meshes (i.e. ten meshes in GDLs, MPLs, CLs, and membranes) are evaluated in figure 3 based on the following equation: where θ ref represents the variables obtained with a strict convergence criterion (i.e. convergence = 10 −9 in the present study). When the convergence criteria are smaller than 10 −3 , the changes in the variables presented in figure 3 are negligible; meanwhile, a further decrease of the convergence criteria will lead to a longer calculation time. Thus, 10 −3 is a desired convergence criterion in the present study as a compromise of modeling accuracy and computational speed. Similarly, a mesh independence study is performed for the ten variables under investigation. The meshes in GDLs, MPLs, CLs, and membranes are divided from 1 to 10 meshes, and the effect of mesh numbers on the numerical results is shown in figure 4. The relative error is calculated based on equation (14), and the reference values are taken from fine-mesh solutions (i.e. ten meshes per layer). As can be seen, finer meshes yield better solutions for all ten variables. The relative errors in temperature, pressure, and concentrations of hydrogen, oxygen, and total species are very small, which is much smaller than 0.1% even for one mesh per layer in comparison with fine-mesh solutions. However, the relative errors in voltage, the dissolved water content in ionomers, and water vapor in porous electrodes are relatively large, and that in liquid water can be as large as 15% for coarse meshes (i.e. one mesh per layer). The findings suggest that if liquid water is not considered, the course-mesh solution can yield comparable results in contrast with the fine-mesh solution, and this coincides with many 1D numerical studies which employed only one mesh per layer [12,33]. It should be pointed out that the mesh number has a significant effect on the liquid water saturation in porous media (see figure 4(j)). For example, the relative error in liquid saturation in anode GDL can be 15% between coarse-and fine-mesh solutions.

Results and data analysis
Two groups of datasets from Toyota are selected to validate the numerical results: low anode relative humidity (RH) (30%) and high anode RH (100%) conditions. The low RH case is operated with anode stoichiometry ratio of 1.5, cathode stoichiometry ratio of 1.25, anode pressure of 2.8 bar, cathode air pressure of 1.1-2.2 bar, and temperature of 358 K; and the high RH case is operated with anode stoichiometry ratio of 1.5, cathode stoichiometry ratio of 1.25, anode pressure of 2.8 bar, cathode air pressure of 1.4-2.5 bar, and temperature of 328 K. It is seen from figure 5(a) that at low anode RH conditions (30%), both single-and two-phase flow models yield good agreement with Toyota data, and the difference in the simulation results between the two models is negligible. This finding is reasonable as the PEM fuel cell is operated under relatively dry conditions, where liquid water does not exist in fuel cell components based on the results obtained from the high-fidelity two-phase flow model. Since the cathode RH is zero and anode RH is 30%, it is difficult for the water vapor to be condensed to liquid water in the fuel cell electrodes at the given steady-state operating conditions. When the anode RH is increased to 100%, the single-phase flow model, in which liquid water is assumed absent, overestimates the voltage under the various current densities from 0.1 to 2 A cm −2 investigated (see figure 5(b)). In contrast, the two-phase-flow model shows good agreement with the Toyota data especially when the current density is larger than 1 A cm −2 , as the liquid water can either block the pores or occupy reaction sites, leading to performance deterioration particularly at high current densities. Figures 5(c) and (d) demonstrate that the ohmic resistance predicted by both two models agrees well with the Toyota data, indicating that the model has been properly validated.
Under low anode RH conditions, the water activity (or RH) is much smaller than 100% (see figure 6). When the PEM fuel cell is operated at a very high current density of 2 A cm −2 , the production water is generated at a relatively high rate due to oxidation-reduction reactions in cathode CLs. In cathode CLs, about 57.11% of the dissolved water generated by electrochemical reactions is transformed into water vapor via phase change, and the accumulation of water vapor increases the amount of water vapor in pores. The cathode channel is fully dried by the inlet flow, leading to a water vapor concentration difference between cathode CL and flow channel. Due to the concentration difference, water vapor is transported toward the flow channel mainly by diffusion. It should be mentioned that at the given conditions, convection is insignificant, and its contribution to the transport of water vapor is negligible. Meanwhile, the remaining 42.89% of dissolved water generated in cathode CLs is transported toward the anode CLs after passing the membranes, and in the anode CLs, dissolved water is transformed to water vapor by phase change and accumulated in the void regions. The local RH (a.k.a. water activity in literature), which is related to water vapor concentration and temperature, becomes slightly larger than 30%. Based on the modeling results, the local RH in anode CL, MPL and GDL is 34.6% 34.3% and 33.7%, respectively. The excess water vapor is transported toward the anode flow channel mainly by diffusion, and no net liquid water production can be found in anode electrodes.
Under a high anode RH of 100%, the water transport is illustrated in figure 7. The cathodic reaction generates a certain amount of water, existing in the ionomer phase of CLs. A portion of dissolved water is transported toward the anode CLs via membranes, and dissolved water in anode CLs is transferred to water vapor. The water vapor accumulated in the pores, leading to a higher local RH (or water activity) which is larger than 100%. Therefore, the excessive moisture is condensed to liquid water, accumulated in electrode pores. When the current density is small (e.g. 0.1 A cm −2 ), about 24.04% of the production water is transported to the anode side with a flux of 0.0125 kmol cm −2 , leading to a slightly higher local RH (>100%) in anode CLs, MPLs, and GDLs. The liquid water transfer rate is correspondingly small, resulting in a low degree of water flooding in the anode electrode. About 5.26% of pore volume in anode CL, 5.23% in anode MPL, and 4.63% in anode GDL are occupied by liquid water. When the current density is high (e.g. 2 A cm −2 ), about 27.99% of the production water is diffused to the anode side with a flux of 0.29 kmol cm −2 . Therefore, the local RH becomes larger, which enhances the liquid water transfer rate. The liquid volume in anode CL, MPL, and GDL becomes 12.9%, 12.8%, and 11.3%, respectively. The water flooding at high current density conditions leads to a severe performance deterioration as shown in figure 5(b). Figure 8 presents the ratio of mass fluxes of hydrogen, oxygen, and water vapor at the interfaces between flow channels and GDLs at different current densities under high anode RH and dry cathode conditions. It is seen that diffusion plays a dominant role in the transport mechanisms of gas species in fuel cells. For instance, the convective mass flux of hydrogen is about 4%-12% of its diffusive compartment (see figure 8(a)), while the ratio of convective and diffusive oxygen fluxes is only 3.5%-5.3% (see figure 8(b)). Under the given conditions, the diffusion dominates at the cathode GDL surface (see figure 8(d)), while at the anode GDL surface, the convective flux can be as large as 20%-35% (see figure 8(c)). Similar results were   reported in Chen et al [58], Sohn et al [59], and Zhao and Li [26]'s work. As can be seen, convection also plays a significant role in the transport of gas species, including hydrogen, oxygen, and water vapor. The transport mechanisms can affect the concentration of fuels and oxidants in CLs at high current densities, which can lead to a substantial increase in activation overpotentials. Meanwhile, the water vapor concentration, which is determined by both convection and diffusion, can lead to different liquid water generation rates in the pores of electrodes, resulting in a significant increase in transport resistances. Most 1D and pseudo-2D models assume that diffusion is the dominant mechanism of mass transport in electrodes by neglecting convection in their governing equations [32,33,37,60]. However, based on the results presented in figure 8, the diffusion-dominant, reduced-dimensional models may not produce a good prediction of fuel cell performance over the wide range of operational conditions experienced in practical fuel cells.
The computational speed of the developed models is shown in figure 9. The program is implemented using Matlab R2016a based on Intel® Xeon® CPU E5-2640 v2 @ 2.00 GHz, RAM 8 GB in the Windows 10 environment. For fine-mesh (ten meshes per layer) domain, the computing time for the two-phase flow model is 1.7 s with a strict convergence criterion of 10 −8 . The computing time becomes 1.1 s when the convergence criterion is reduced to 10 −3 (a desired criterion based on modeling accuracy as discussed early), indicating about a 35% drop of computing time. If liquid water is not computed (i.e. the single-phase flow model), the computing time is reduced to 0.69 s-0.48 s for the convergence criteria of 10 −8 and 10 −3 , respectively. For coarse-mesh (one mesh per layer) domain, the computing speed is much faster. For example, the computing time required for the steady-state single-phase flow model is only 0.26 s, while that of the two-phase flow model is about 0.56 s.

Conclusions
In this study, computationally efficient PEM fuel cell performance models considering single-and two-phase flows are developed based on a SIMPLE algorithm with a hybrid central and UDS. The models take various transport phenomena and electrochemical behaviors into account, including mass conservation, momentum transport, species transport, liquid water, membrane water, Nernst equation, and simplified Butler-Volmer equation. Both single-and two-phase flow models exhibit excellent agreement with the experimental data in terms of polarization curve and ohmic resistance for low relative humidity conditions. For high relative humidity operation, the single-phase-flow model over-estimates the output voltage of PEM fuel cells under high current density conditions as liquid water is not considered. When liquid water is incorporated in the two-phase flow model, the modeling results agree well with the experimental data. The modeling results also suggest that the maximum convective/diffusive ratio of H 2 , O 2 , and vapor mass fluxes can be 12%, 5.3%, and 35%, respectively, which are ignored in most diffusion-dominant models. The ignorance of convection mechanism in gas species transport within porous electrodes may lead to inaccurate modeling results if the electrode diffusivity is measured experimentally. In addition, both models show excellent computing speed. For example, for the coarse-mesh domain, it requires only 0.26 s and 0.56 s for the single and two-phase flow models to converge with the desired convergence criterion of 10 −3 . These computationally efficient models will be beneficial for the control of the fuel cells in actual automotive and industrial applications. It should be mentioned that the 1D model still has certain limitations due to the nature of 1D models. The variations in reactant concentration and temperature in channels are commonly assumed negligible in 1D models. When the fuel cell stack is very large, the non-uniform distribution of reactants and temperature may introduce errors in the boundary conditions of 1D models, which should either be corrected or use 1 + 1D models (or pseudo-2D models) to address these variations. Further, it is important to explore the transient high-fidelity 1D models to evaluate the dynamic behaviors for real-time control applications.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).