Free open reference implementation of a two-phase PEM fuel cell model

In almost 30 years of PEM fuel cell modeling, countless numerical models have been developed in science and industrial applications, almost none of which have been fully disclosed to the public. There is a large need for standardization and establishing a common ground not only in experimental characterization of fuel cells, but also in the development of simulation codes, to prevent each research group from having to start anew from scratch. Here, we publish the first open standalone implementation of a full-blown, steady-state, non-isothermal two-phase model for low-temperature PEM fuel cells. It is based on macro-homogeneous modeling approaches and implements the most essential through-plane transport processes in a five-layer MEA. The focus is on code simplicity and compactness with only a few hundred lines of clearly readable code, providing a starting point for more complex model development. The model is implemented as a standalone MATLAB function, based on MATLAB's standard boundary value problem solver. The default simulation setup reflects wide-spread commercially available MEA materials. Operating conditions recommended for automotive applications by the European Commission are used to establish new fuel cell simulation base data, making our program a valuable candidate for model comparison, validation and benchmarking.


Introduction
The development of macro-homogeneous models of the membrane electrode assembly (MEA) of low-temperature proton exchange membrane fuel cells (LT-PEMFCs) goes back almost 30 years, to Springer et al. [1] and Bernardi & Verbrugge [2][3][4]. The first non-isothermal variant was published by Fuller & Newman [5]. Ever since these pioneering efforts, research and development on numerical simulations of the various transport processes within the MEA have brought forward a large variety of fuel cell models at different length scales, helping scientists and engineers to better understand the complex nonlinear behavior of these promising energy converters.
Even though many of those models are based on the same core functionality, the policy of publishing the mathematical model description but keeping the numerical implementation tightly closed, has forced software developers to reinvent the wheel by starting from scratch over and over again. The unavailability of a fully transparent and easy-to-understand reference implementation of a basic MEA model with spatial resolution is slowing down the advent of modeling in the fuel cell community, some participants of which even hesitate to include modeling in their work altogether. In their comprehensive review article [6], Weber et al. note that "the majority of PEFC models to date have either been implemented in commercial software such as FLUENT, COMSOL, STAR CD, or implemented in-house. In either case, the source code has not been made available to the public. This has several major drawbacks including (i) lack of validation and comparison between models, (ii) lack of extension capabilities, and (iii) implementation limitations". * Corresponding Author: roman.vetter@zhaw.ch They also comment that "the key disadvantages of most open-source codes are no graphical user interface and a necessary knowledge of Linux OS and [...] C++ or python".
Open-source code development and validation activities for fuel cells are only recently picking up some steam. The International Energy Agency has launched an annex on open-source modeling of fuel cells systems, but the focus has primarily been on solid oxide fuel cells so far [7,8]. To this day, there are only two known open-source codes capable of simulating the state of the art in PEMFC modeling at the cell scale: • OpenFCST [9], a C++ package based on the deal.II finite element library, freely available under the MIT license. It is highly capable, but with more than 120 000 lines of C++ code as of version 0.3, it is also very heavy-weight and difficult to handle. • FAST-FC [10], a finite volume tool built on top of Open-FOAM. It consists of about 12 000 lines of code (not counting OpenFOAM) that are published under the GNU General Public License v3, which can pose an insurmountable legal barrier for commercial use. A third open-source toolkit that has been used to study porous media of PEMFCs is OpenPNM [11], a pore-network model implementation in Python/SciPy. It is freely available under the MIT license and consists of about 25 000 lines of code.
With this paper, we present a very light-weight, free standalone implementation of a full-blown macro-homogeneous five-layer MEA model for low-temperature PEM fuel cells. The model is formulated in three dimensions, but implemented only in one spatial dimension to represent the dominating through-plane transport processes. It is nonisothermal and two-phase to capture important thermal effects such as phase-change induced flow, but isobaric and stationary to avoid the computational complexity arising from pressure gradients and unsteady behavior. We designed the program code for • Simplicity and compactness. It comprises less than 400 lines of commented code and does not require manual differentiation of the equations. This dramatically simplifies modifications such as the substitution of individual parameterizations or boundary conditions, which can be done by replacing a single line of code. No lookup tables or data interpolation are involved. • Portability and compatibility. The model is implemented as a standalone MATLAB function, relying only on MAT-LAB's standard boundary value problem solver. This choice of programming environment lets the simulation run on a large variety of platforms and also lets it benefit from MATLAB's widespread availability in science and industry. • Transparency. The model equations and boundary conditions are fully disclosed, and the complete simulation output is shown in the paper, including all potentials, fluxes and the entire polarization curve. • Accessibility. The program is well documented and ready to be used out of the box. Modifying the code requires only minimal programming knowledge, running it requires none at all. All plots are automatically generated.

Mathematical model
In the one-dimensional PEMFC model developed here, the MEA is represented as a series of five adjacent homoge-neous interval subdomains representing a cell-area-averaged through-plane section of the porous layers of a fuel cell. It includes the proton exchange membrane (PEM) in the middle, sandwiched by two catalyst layers (CLs) and two gas diffusion layers (GDLs). The gas channels (GCs) and monoor bipolar plates on either end are modeled as boundaries of the MEA. Other subdomains such as microporous layers (MPLs) are not explicitly modeled, but can be added without difficulty. The geometry of the MEA model is shown in Fig. 1. We keep the mathematical description as compact as possible, briefly summarizing the conservation laws and transport equations of the model.

Electrochemistry.
In hydrogen-fueled LT-PEMFCs, the net electrochemical reaction is and hence the reversible cell potential is given by the Nernst equation [13] where F is the Faraday constant, R the gas constant, P ref = 1 atm is the reference pressure, T is the temperature, pH 2 = xH 2 P and pO 2 = xO 2 P are the partial pressures of hydrogen and oxygen and ∆G = ∆H − T ∆S is the Gibbs free energy change of the reaction. We now assume that the overall redox reaction can be split into a single-step hydrogen oxidation reaction (HOR) in the ACL and a single-step oxygen reduction reaction (ORR) in the CCL, and that both can be described with sufficient accuracy by Butler-Volmer kinetics. The reaction rate in the homogenized catalyst layers is thus locally given by [13][14][15] with (positive) activation overpotential ∆φ = φe − φp is the Galvani potential difference between the electron and proton conducting phases (see Tab. 1), and (5) with total reaction entropy ∆S = ∆SHOR + ∆SORR. The sign convention used here is such that a positive i corresponds to a source of positive charge or mass in the continuity equations (see Tab. 2).
Transport of charge, heat and mass. For the remainder of the model description we follow the continuum approach to describe the most dominant transport processes of charge, energy, gas species and water by conservation laws. This results in eight coupled second-order partial differential equations (PDEs), which are summarized in Tab. 1 for the steady state.
In the CLs and GDLs, Ohm's law is assumed to govern how the flux of electrons je is driven by a gradient of the electronic phase potential. The analogous equation is used for the flux of protons jp within the electrolyte phase of the CLs and the membrane. The two electrostatic phase potentials φe and φp coexist in the CLs (see Fig. 1), defining ∆φ in these two domains. The approach is adopted from the classical porous-electrode theory of Newman, where it is assumed that the electric double layer constitutes only a small volume compared to any of the phases or the electrode itself [16].
Heat conduction is the dominating mode of energy transport in the MEA [17], allowing for an accurate description of the heat flux jT by Fourier's law in all five subdomains. This is the third differential equation.
The description of water balance in the ionomer is based on the seminal model by Springer et al. [1]. To represent the degree of humidification, the number of water molecules per acidic group λ is used. The molar flux of dissolved water j λ is composed of the sum of its two most significant contributions: back diffusion due to a moisture gradient (j λ ∼ −∇λ) and electro-osmotic drag (j λ ∼ jp ∼ −∇φp).
The next three equations are dedicated to the transport of gas species on both sides of the membrane. If gas crossover is neglected, it is sufficient to consider hydrogen only on the anode side and oxygen on the cathode side, whereas water vapor is present in both gas mixtures. A third gas component is implicitly accounted for on either side (typically nitrogen if air is supplied to the cathode) and need not explicitly be computed because the sum of mole fractions X xX = 1 everywhere. We assume uniform gas pressure P in the steady state, and hence the dominant transport mechanism is inter-diffusion of the gas species. The simplest transport model for this is Fick's law jX = −CDX ∇xX [18], which is employed here for all species. Thermal diffusion, as it results from the chemical potential gradient as the general driving force for species transport, is neglected. The ideal gas law is used to calculate the interstitial gas concentration C = P/RT .
Finally, for the description of liquid water transport, we adopt the unsaturated flow theory that was carried over from soil physics to the fuel cell modeling community by Natarajan & Nguyen [19], and which has since become the de-facto standard in macro-homogeneous two-phase MEA modeling. Darcy's law is transformed into an equation for liquid water flux driven by a gradient ∇s, where s denotes the liquid water saturation (fraction of pore space filled with liquid water). This requires the specification of both the saturation-dependent hydraulic permeability κ and the differential relationship between the capillary pressure and saturation, ∂pc/∂s, as material properties.
For each of these eight fluxes, a continuity equation is expressed in the last column of Tab. 1, equating the divergence of each flux j with a corresponding source term S. These eight PDEs become nonlinear when the coefficients and/or source terms are expressed in terms of the dependent variables. All phase transitions and reaction rates appear as sources that couple these PDEs as detailed in the following.
Source terms and phase transitions. A summary of all source term definitions is given in Tab. 2. In the ACL, hydrogen is split into electrons and protons with a reaction rate given by Eq. 3, giving rise to source terms Se and Sp. Faraday's law determines the rate of hydrogen consumption in the ACL as well as the oxygen consumption and water production in the CCL (moles consumed or produced per unit volume, time and exchanged electron pair): The molar oxygen consumption rate is only half this much (see Eq. 1). It is assumed that water is produced at the platinum-ionomer phase boundary in dissolved form [20], hence SF appears as a contribution to S λ in the CCL. Absorption and desorption of water vapor into/from the ionomer does not happen instantaneously, but at a finite rate over a time span of the order of an hour [21][22][23][24]. To account for this significant ionomer-gas interfacial water transport resistance, the sorption source term S ad appearing in the continuity equations of λ and xH 2 O is set to [25] S ad = where L is the thickness of the CL, Vm the molar charge volume of the ionomer, λeq denotes the RH-dependent equi- librium water content of the ionomer, and ka, k d are materialdependent mass-transfer coefficients. The commonly employed approach to include liquid-vapor phase change in macro-homogeneous MEA modeling is to assume the mass transfer to be driven by the vapor partial pressure difference to the saturation pressure [26]. With xsat = Psat/P , the corresponding water source/sink term can be expressed as where γe and γc are the evaporation and condensation rates. The latent heat released or absorbed during the above two phase transitions can be modeled by adding the following contributions to the total heat source ST : H ad and Hec are the molar enthalpies of desorption and vaporization. Joule's first law provides two more heat sources induced by the electric and ionic currents, And finally, the heat dissipated by the electrochemical reaction is given by the sum of activation and Peltier heats [6]: Boundary conditions. To complete the mathematical model description, a set of boundary conditions (BCs) needs to be specified (two for each contiguous support in each of the eight second-order PDEs), which are listed in Tab. 3. The membrane is assumed to be impermeable for all gas species as well as for electrons and liquid water, hence these normal fluxes vanish at the membrane boundaries. Protons and dissolved water on the other hand are bound to the ionomer phase, which implies zero fluxes at the outer surfaces of the catalyst layers. Since the electrostatic potentials can be freely offset, φe is set to zero at the outer AGDL boundary. At the remaining interior subdomain interfaces, continuity of the potentials and fluxes is assumed as indicated in Tab. 3. For the remaining outer boundaries, a reasonable choice depends on the scenario to be simulated. We impose the cell voltage U by applying the Dirichlet BC φe = U at the other end of the MEA, but equivalently one can control the current density I by using a Neumann BC n · je = I instead. We furthermore impose the gas channel temperature, pressure, relative humidity and gas composition via where αH 2 (αO 2 ) is the hydrogen (oxygen) mole fraction in the supplied gas when dry. Perhaps the most delicate interface treatment is that of liquid water, and the formulation of physically accurate models is a topic of ongoing research in PEMFC modeling [6]. At the CGDL/GC interface, water droplets form and detach in a dynamic fashion (e.g., [27,28]) that is difficult to translate into a steady-state area-averaged BC. Historically, simple Dirichlet BCs for s are often used [19,29,30]. Conditional unidirectional flux conditions have later been proposed as a more realistic replacement [31,32], but solving these numerically can be a challenge [33]. We use here a Dirichlet BC for s at the CGDL/GC interface-the simplest common denominator in two-phase MEA modeling-bearing its limitations in mind.
Initial conditions. Nonlinear problems require a good initial guess of the solution for iterative solvers to converge. It is most convenient to iterate over cell voltages from high to low to generate the polarization curve, since one can then start with all-zero fluxes as a good initial guess. For the potentials, the following initial conditions are usually sufficient for convergence:

Parameterization
For the establishment of a useful reference PEMFC simulation suitable for model comparison and benchmarking, it is important to furnish the model with well-established parameterizations of widely available commercial MEA materials. Nafion NR-211 is the membrane of choice here, owing to the large market share and the vast pool of characterization data of Nafion in the literature [34]. Since Toray carbon paper is among the most comprehensively characterized GDLs in the literature, we use Toray TGP-H-060 material properties to parameterize the GDLs in the model. Together, these materials form a typical modern MEA as it may, for instance, be used for automotive applications. Standard literature data are used for the remaining material-independent electrochemical and physical properties.
Water properties. For water produced in liquid form at 25 • C and 1 bar, the enthalpy of formation is ∆H = −285.83 kJ/mol [35]. The saturation pressure of water vapor Psat can be approximated with the Antoine equation in the temperature range T = 50−100 • C [36]: The same relationship goes by the name of Vogel equation when used for the dynamic viscosity of liquid water µ, and it has the following coefficients in the temperature range The condensation and evaporation rates may be estimated as [20,29] γe = kea lg s red where a lg ≈ 2 m 2 /cm 3 is an effective liquid-gas interfacial surface area density scaling factor [20], and ke and kc are the Hertz-Knudsen mass transfer coefficients, given at atmospheric pressure by [38] ke where Mw = 18 g/mol is the molar mass of water. In Eq. 15, liquid water saturation dependence of the phase change interface is introduced through the reduced saturation where sim denotes the immobile or inaccessible saturation (i.e., liquid water that does not contribute to transport pathways or phase change, e.g. due to spatial isolation). It is estimated as sim = sC. Finally, the latent heat of evaporation/condensation is Hec ≈ 42 kJ/mol in the temperature range relevant for PEMFC operation [39].
Electrochemical parameters. Since the ORR is the ratelimiting half-reaction, it is crucial to model the concentration and temperature dependence of the exchange current density in the cathode with high accuracy. Neyerlin et al. [40] have obtained the following relationship for a Pt/C cathode: Albeit measured for a low equivalent weight (EW) ionomer, the activation energy in Eq. 18 is consistent with reported values for higher EW ionomers (such as 1100 EW Nafion) as well [40]. Furthermore we assume that Eq. 18, which was fitted premising the Tafel equation, can be applied to the Butler-Volmer equation without modification. For the much faster HOR, the Butler-Volmer equation has been reported to hold with [41,42] We consider a cathode platinum loading that is three times as high as in the anode: a = 1×10 11 cm 2 Pt /m 3 in the ACL and a = 3×10 11 cm 2 Pt /m 3 in the CCL. The symmetry factor β is assumed to be 1/2 in both half-reactions. Lampinen & Fomino's values for the half-reaction entropies ∆SHOR = 0.104 J/mol K and ∆SORR = −163.3 J/mol K [43] appear to be the most plausible in the literature and are therefore adopted here.
Ionomer-related parameters. For Nafion, the enthalpy of (de-)sorption is almost equal to that of vaporization when not completely dry [54], hence we set H ad = Hec. For the ionic conductivity of Nafion membranes, a power law from percolation theory with Arrhenius temperature correction was found to best fit various experimental data [55]: where denotes the volume fraction of water in the ionomer. Vm = 1020/1.97 cm 3 /mol is the equivalent volume of the dry membrane (EW [56] divided by mass density [45]) and Vw = 18/0.978 cm 3 /mol the molar volume of liquid water at typical PEMFC operating conditions. In Eq. 20 the Bruggeman correction 1.5 i is used to account for the different ionomer contents i in the PEM and CLs [57].
Another crucial transport parameter is the water diffusivity in the ionomer D λ . Experimental difficulties have led to large discrepancy in the literature data for Nafion [34]. The measurements carried out by Mittelsteadt & Staser [58] appear to be among the most sophisticated. We refitted their data for Nafion membranes by a rational polynomial in λ to obtain a smooth parameterization that captures all  [53]. essential features of the data: This new fit for D λ is plotted in Fig. 2. Analogous to Eq. 20, the Bruggeman correction is used to model water diffusion through the partial ionomer content of the CLs. For the electro-osmotic drag coefficient, we adopt Springer's original linear law [1] and their sorption isotherm at T = 30 • C is used to determine the equilibrium water content of the ionomer in Eq. 7: where RH = xH 2 O/xsat is the relative gas humidity. A parameterization for the mass-transfer coefficients of water vapor in Nafion membranes has been determined by Ge et al. [59]: where aa = 3.53×10 −3 cm/s and a d = 1.42×10 −2 cm/s. Transport in the porous media. To estimate the Fickean gas diffusivities within the pore space of the partially flooded CLs and GDLs, we amend the Chapman-Enskog formula [18] by the usual porous media correction factor for porosity and tortuosity, as well as by a saturation correction (1 − s) raised to the third power [60]. For X = H2, O2, H2O, with the following prefactors: • As for the capillary pressure-saturation relationship, a large number of correlations have been suggested for various GDLs [61]. We employ here the following expression that has been determined for Toray TGP-H-060 specifically [62]: The second parameter in the liquid water flux equation is the hydraulic permeability κ, which has a great impact on the liquid water distribution. It is modeled as [29] κ = 10 −6 + s 3 red κ abs (28) where κ abs denotes the absolute (intrinsic) permeability of the porous medium and the small numerical tolerance is added to avoid the singularity at s red = 0. A reasonable estimate for the liquid water saturation at the CGDL/GC interface sC can be found with the notion of equivalent capillaries in the GDL through which the liquid water is transported. Using the Young-Laplace equation pc = 2γ cos θ/r with effective surface contact angle θ = 130 • [62] and equivalent capillary radius r = 40 µm [63] for Toray TGP-H-060, one finds an equivalent capillary surface pressure of pc ≈ 2 kPa. Using Eq. 27, this translates to sC ≈ 0.12, which is the boundary value used here. Note that this does not pose any limitation on the boundary flux js at the CGDL/GC interface, i.e., the interfacial liquid water flux will automatically be such that this pressure BC is met.
All remaining material properties are considered constant as listed in Tab. 4. The GDLs are assumed to be moderately compressed from 190 to 160 µm, corresponding to an applied clamping pressure of about 1.4 MPa [44].

Numerical implementation
The model is implemented as a standalone MATLAB function MMM1D (short for one-dimensional Master MEA Model)  The complete source code as printed in the Appendix can be obtained from https://www.isomorph.ch for free. It is released under the 3-clause BSD license available from the same website, permitting unrestricted commercial and noncommercial use subject to the condition that the original source be referenced.
In the provided reference implementation, moderate bvp4c error tolerances are used (relative: 10 −4 , absolute: 10 −6 ), resulting in an average of 54 mesh nodes used in total. The execution time for a full sweep over all cell voltages in steps of 50 mV is a few seconds on a modern laptop computer, and the maximum absolute (relative) discretization error under these conditions is 4.4 mA/cm 2 (0.23%). Increased accuracy can be obtained by reducing these tolerance values if desired. All output plots shown in this paper have been obtained with very high accuracy using lower error tolerances (relative: 10 −6 , absolute: 10 −10 ), which resulted in 158 mesh nodes on average for the base case.

Simulation results
Base case. By default, our program code simulates a PEMFC at typical operating conditions referred to as the base case. These conditions are listed in Tab. 5 and are used here to present the complete simulation output.
The polarization curve is plotted in Fig. 3 for the entire range of cell voltages in steps of 10 mV. Fig. 4 shows the eight potentials and Fig. 5 the corresponding fluxes across the MEA layers as predicted by the model at different cell voltages in steps of 100 mV. Each subplot is restricted to the support of the respective variable for a more detailed view where possible. All these plots are automatically generated by the MATLAB function.
The membrane phase potential and water content profiles deserve closer attention. Even though Sp ≡ 0 in the bulk membrane, φp(x) exhibits significant curvature (Fig. 4). This is due to a strong spatial variation of the proton conductivity through the membrane, caused by a relatively steep decline of λ toward the anode, which in turn is the result of strong electro-osmotic drag of dissolved water to the cathode. Parameterizations of the water diffusivity that predict larger values than Eq. 22, lower drag coefficients than Eq. 23, and higher ionic conductivities than Eq. 20 at low water content all yield higher λ near the anode. This results in more flat potential profiles φp(x) and consequently, higher current densities. These parameters are, in fact, among the material properties with the largest impact on the predicted fuel cell performance [65]. An extensive study on this subject is currently underway at our institute.
To further extend the data basis for future model comparison, we report on a few additional model characteristics for the base case, derived from the model output shown above. Tab. 6 lists some key figures, among which are the ohmic membrane resistance (excluding the contribution from the ionomer in the CLs) the average MEA temperature and the mean water content of the ionomer Here, the integration runs over the whole catalyst-coated membrane (CCM = ACL ∪ PEM ∪ CCL). These quantities are evaluated at U = 0.6 V, but it is straightforward to use our program to calculate them at any other operating point.     EU harmonized stress tests. A common limitation of fuel cell research in science and industrial applications to date is the lack of international standardization for experimental characterization and numerical modeling. To this end the Joint Research Centre (JRC) of the European Commission has recently issued a set of single-cell stress tests [12], termed T1 to T9 here, of which the first seven are directly applicable to the present 1D model. They define operating points aimed at characterizing the performance of a PEMFC under automotive conditions with variations in temperature (T1 and T2), gas humidification (T3 to T5) and gas pressure (T6 and T7), as detailed in Tab. 7. In order to establish a new baseline for harmonized fuel cell model comparison, validation and benchmarking, we subject our model to this series of stress tests and report the normalized output in the form recommended by the JRC in this section. These results may serve for comparison purposes in future PEMFC model development efforts.
The performance criteria agreed upon by the JRC specify three points on the polarization curve to be evaluated for each stress test. Tab. 7 lists the model output data for the reference case and the applicable stress tests. The full resulting polarization curves are juxtaposed in Fig. 6. As proposed by the JRC, the three measured values are normalized by the corresponding values at reference conditions according to normalized result = 1 − reference result stress test result (32) and visualized as a spider plot in Fig. 7. Since the limit current density of T2 is only 0.733 A/cm 2 , the cell voltage of T2 at 0.8 A/cm 2 is omitted from this compilation. Only Tests T1-T7 were conducted, because T8 and T9 (stoichiometry variation) are inapplicable to the present 1D model. When interpreting the stress test results, it is important to note that a 1D through-plane MEA model simulates the characteristics of a differential fuel cell, i.e., one with a very small active cell area. The JRC specifies gas inlet conditions irrespective of gas channel length or cell area. In fuel cells with sufficiently long channels, the supplied gas gets humidified significantly as it flows downstream in the flow channels, making the area-averaged relative humidity much larger than at the inlet. With RH values between 20 and 85%, the JRC test specifications therefore represent very dry conditions when applied to a differential cell. In order to mimic typical area-averaged conditions in the base case, the relative humidities were set to larger values for the default simulation setup (see Tab. 5), whereas we strictly abide by the dry JRC specifications here. Only T1 is humid enough to allow for the presence of liquid water in a differential cell as predicted by the model. All other stress tests are too dry for the RH to reach 100% within the MEA of a differential cell. We therefore set sC = 0 for all stress tests but T1.
The performance results shown in Fig. 6 are thus strongly moisture-limited. T1, being the most humid scenario, allows for the best performance at moderate current densities, before electro-osmotic drag dries out the anode end of the membrane and the low temperature of only 45 • C becomes the limiting factor in the ionic conductivity (cf. Eq. 20), allowing the drier but hotter T3-T5 to perform better at low cell voltages. T2 ("dry/hot desert") yields the worst performance prediction due to the extremely dry membrane (λ ≈ 2.5 on average). Owing to the stronger back diffusion of dissolved water in the ionomer, T4 (dry fuel) slightly outperforms T3 (dry air) in spite of the better ACL humidification under T3. Finally, we note that the dry test conditions of T6 and T7 do not allow the differential fuel cell to reach the large current density regime where a lot of fuel or oxidant is consumed and starvation becomes an issue. Under these test conditions, the model therefore exhibits very little sensitivity to pressure variations. This characteristic changes when less permeable (thicker) diffusion media and more humid gases are employed, which can easily be verified by using the provided MATLAB program.

Conclusion
With our free open MATLAB implementation of a macrohomogeneous two-phase PEMFC model, anyone with access to a MATLAB installation can readily run a state-of-the-art PEMFC simulation, substitute material parameterizations, add new model features or conduct parameter studies. It offers several major advantages over existing open-source codes. With less than 400 lines of compact, commented MATLAB code, it is exceptionally easy to read, maintain, extend and embed in other programs, with no third-party software, compilation of source code or knowledge of Linux and C++ or Python required. Both commercial and noncommercial use are permitted thanks to the BSD-like licence under which the program is published. We have provided the complete simulation output of the model at typical operating conditions, laying a new cornerstone in the ongoing effort of making numerical PEMFC modeling more transparent and accessible.
Despite these evident strengths, our model-like every model developed thus far-has limitations that need to be kept in mind. For the sake of simplicity, we intentionally neglect several chemical and physical processes occurring in real fuel cells that might become relevant under certain operating conditions. Among others certainly, the model does not account for the effects of heat convection, gas convection (due to a non-zero gas pressure gradient), gas species permeation and pressure-driven hydraulic permeation of water through the membrane, thermo-osmosis, Schroeder's paradox [66], Knudsen diffusion, electrical and thermal contact resistance [67], mechanical deformation and other effects of clamping pressure, non-uniformity in material properties such as wettability and porosity, double layer effects, multistep reaction kinetics, platinum oxide formation and the change of Tafel slope [6], non-uniformity of ionic concentrations, ionic species migration, water droplet formation and detachment at the GDL/GC interface, the short-range effect observed in thin porous layers [68,69], degradation, unsteady phenomena, gravity, ice formation and melting, etc.
Many of the parameters and transport coefficients adopted for the present MEA model are subject to relatively wide variation in the literature and between different materials, while having a significant impact on the simulation results at the same time. Perhaps the largest source of uncertainty in the model lies in the liquid water flux through the porous media and across the interfaces. These are indeed modeling aspects for which no satisfactory universal solutions exist to date [6]. While PEMFC model development is still an on-going process, an accessible numerical tool with which new parameterizations, interface conditions etc. can easily and quickly be tested, compared or validated against measurement data, can be key to further progress in this direction. Our open reference implementation of a 1D MEA model meets these requirements. A demonstration of how it can be utilized to quickly assess different material parameterizations has recently been presented [65]. With a runtime of about a second on an ordinary laptop computer for a single simulation, it is suited even for time-critical applications. In cases where even less resources are available, it is possible to simplify the model in a number of ways, such as omitting the explicit account for liquid water through artificial extrapolation of Eq. 24 to the supersaturation regime as in ref. [1], merging the two diffusion equations for hydrogen and oxygen into one (because they are solved on disjunct subdomains), or removing the gas transport equations altogether. Moreover, as the first plot in Fig. 4 shows, the electron phase potential φe varies only little through the cell depth in the simulated base case. An order-of-magnitude analysis shows that for a GDL with thickness L ∼ O(100 µm) and electric conductivity σe ∼ O(10 3 S/m), the voltage loss associated with a current density of I ∼ O(1 A/m 2 ) going through it is IL/σe ∼ O(1 mV). In cases where voltage drops in this order of magnitude and the corresponding ohmic losses of I 2 L/σe ∼ O(1 mW/cm 2 ) are deemed insignificant, the potential φe may be replaced by constants in the GDLs. Doing so also in the CLs, on the other hand, would violate charge conservation (Se = −Sp, cf. Tab. 2) and lead to convergence difficulties.
Conversely, there is also much room for model extensions. Aside from the inclusion of the above-mentioned neglected effects, the model can be augmented by adding additional subdomains to represent MPLs, by using the Brinkman equation in place of Darcy's law, the Maxwell-Stefan equations in place of Fick's law for gas diffusion, the Nernst-Planck equation in place of Ohm's law, etc. Moreover, it is straightforward to add liquid water also on the anode side if desired. It is also possible to deeply refine the model in terms of its parameterization, for instance by including temperature dependence in the water sorption isotherm of the ionomer or in the electro-osmotic drag coefficient. Additional model complexity, detailed material property parameterizations and higher-dimensional models are being developed at our institute and are available upon request. More information can be found at https://www.isomorph.ch.