{'fraction_non_alphanumeric': 0.03868929469366151, 'fraction_numerical': 0.020725817123197036, 'mean_word_length': 4.687394167450612, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Experimental parameter uncertainty in PEM fuel cell modeling. Part I: Scatter in material parameterization

Ever since modeling has become a mature part of proton exchange membrane fuel cell (PEMFC) research and development, it has been plagued by significant uncertainty lying in the detailed knowledge of material properties required. Experimental data published on several transport coefficients are scattered over orders of magnitude, even for the most extensively studied materials such as Nafion membranes, for instance. For PEMFC performance models to become predictive, high-quality input data is essential. In this bipartite paper series, we determine the most critical transport parameters for which accurate experimental characterization is required in order to enable performance prediction with sufficient confidence from small to large current densities. In the first part, a macro-homogeneous two-phase membrane-electrode assembly model is furnished with a comprehensive set of material parameterizations from the experimental and modeling literature. The computational model is applied to demonstrate the large spread in performance prediction resulting from experimentally measured or validated material parameterizations alone. The result of this is a ranking list of material properties, sorted by induced spread in the fuel cell performance curve. The three most influential parameters in this list stem from membrane properties: The Fickean diffusivity of dissolved water, the protonic conductivity and the electro-osmotic drag coefficient.


Introduction
Fuel cell researchers that consider numerical modeling for their work are confronted with the difficult question which material properties to measure and plug into the models in order to obtain reliable simulation results. With finite financial and time budgets at hand, effort put into the characterization of membrane-electrode assembly (MEA) components must be prioritized. Is it worthwhile to determine the tortuosity of gas diffusion pathways through the porous layers with high resolution imaging, or should the focus be on precise measurements of water transport through the proton exchange membrane (PEM)? How will uncertainties in these parameters affect the predictive power of a detailed fuel cell model? MEA engineers face similar problems: For the efficient development of improved multifunctional materials for PEM fuel cells (PEMFCs), it is critical to know what quantity to optimize for. Will a thinner membrane yield significant improvements in fuel cell performance, or is it more productive to look into reducing contact resistances? Even when the MEA setup is fixed and known in detail, modelers are confronted with the problem of picking the right material parameterization among the many that have been proposed in the open literature. Published data on the most critical transport processes occurring in PEM fuel cells (NB for the very same materials) are sometimes scattered significantly, which raises the question which to adopt.
Recently, we have demonstrated how difficult it can be to reliably predict fuel cell performance with numerical modeling under the present spread of experimental material data in the literature [1]. The water diffusivity in the membrane [2,3], electro-osmotic drag [4], protonic conductivity of the membrane [5] and the phase change rates of water [6] are some out of several properties that are required for state-of-the-art twophase PEMFC models and for which measured data disagree largely; further ones are discussed in this paper.
Not many studies systematically address the variation in performance prediction of fuel cell models arising from uncertainty in their parameterization. Typically, a single specific constitutive parameterization is implemented for each material property, consisting of a functional relationship and a set of fit parameters. Only the latter have been varied to estimate the model sensitivity to them (e.g., [7][8][9]). Uncertainty resulting from the choice of the complete constitutive relationships, on the other hand, is largely unexplored. In Part I of this bipartite paper series, we intend to fill this gap. After a brief model summary in Sec. 2, we discuss its constitutive parameterization in detail in Sec. 3, reviewing the literature on several essential material properties. With this comprehensive database of proposed material parameterizations at hand, the scatter in fuel cell performance resulting from it will be presented in Sec. 4. This allows us to conclude with a sorted list of MEA properties which cause uncertainty in the resulting polarization curve and for which more reliable experimental characterization is therefore required. In Part II [10], we will extend this survey by carrying out extensive forward uncertainty propagation analyses to get a more rigorous insight into the relative importance of PEMFC material parameterizations.

Fuel cell model
In order to get a quantitative picture of the uncertainty in the performance prediction of state-of-the-art PEMFC models due to disagreement or scatter in the experimental literature data on fuel cell materials and processes, a model is needed that is both easily modifiable and numerically efficient. Here, we build upon our previously developed macro-homogeneous, steady-state, two-phase model of a five-layer MEA [11]. Being built for easy parameter substitution and fast execution of many different simulations (e.g., with different constitutive material properties), this one-dimensional model offers the best middle ground between model complexity, extensibility and computational efficiency. It models the coupled charge, heat and mass transport processes in through-plane direction within a differential PEMFC. For improved accuracy, we extend it here in a number of ways: Among other amendments (detailed further below), multicomponent gas diffusion is accounted for with the Maxwell-Stefan equations, Knudsen diffusion is included, thermo-osmosis and Schroeder's paradox are incorporated, and electrical as well as thermal contact resistance between the MEA layers are added to the model. In this section we briefly recapitulate the governing equations and refer to [11] for additional information.

Conservation laws
The MEA modeled here consist of a classical symmetrical arrangement of five porous layers: A polymer electrolyte membrane (PEM, Nafion) in the middle, sandwiched by two thin catalyst layers (ACL & CCL) containing an electron-conducting carbon matrix and a proton-conducting ionomer phase (Nafion), which in turn are clamped between two fibrous gas diffusion layers made of hydrophobized carbon paper (AGDL & CGDL).
For the conservation of electrons (subscript e) and protons (p), Ohm's law je = −σe∇ϕe, div je = Se is solved with effective layer-dependent conductivities σe and σp (to be specified in Sec. 3.2). The source terms read where the anode and cathode reaction rates SA and SC will be specified in Sec. 3.1. In the CLs, the two phase potentials ϕe and ϕp coexist, defining the total Galvani potential difference ∆ϕ = ϕe − ϕp. Using the reversible half-cell potentials the anode and cathode activation overpotentials are defined as ηA = ∆ϕ − ∆ϕ 0 A and ηC = ∆ϕ 0 C − ∆ϕ.
∆H ref is the standard enthalpy of formation of liquid water, ∆S ref A and ∆S ref C the half-reaction entropies, F Faraday's constant, R the gas constant, T the absolute temperature, P ref the reference pressure, pH 2 = yH 2 P the partial pressure of hydrogen, and pO 2 = yO 2 P that of oxygen.
For the conservation of energy, heat conduction is considered according to Fourier's law across all layers: with effective thermal conductivity k (Sec. 3.3). The total heat source ST includes the following contributions: Ohmic heating in the electron and proton-conducting phases, Peltier heating and thermal activation losses, and the latent heat of water ab-/desorption (subscript ad) in the CLs and evaporation/condensation (ec) where liquid water is present (assumed only on the cathode side of the membrane). This yields j 2 e /σe in AGDL j 2 e /σe + j 2 p /σp + (ΠA + ηA)SA + H ad S ad in ACL j 2 p /σp in PEM j 2 e /σe + j 2 p /σp + (ΠC + ηC)SC + H ad S ad + HecSec in CCL j 2 e /σe + HecSec in CGDL (6) where ΠA,C = −∆S ref A,C T /2F are the Peltier coefficients of the half-reactions, S ad (Sec) the reaction rate of water ab-/desorption (evaporation/condensation) detailed in Sec. 2.2, and H ad (Hec) the corresponding latent heat (Sec. 3.12).
The transport of dissolved water through the membrane and the ionomer phase of the CLs is modeled in the spirit of Springer et al. [12], who simplified the general diffusion equation with chemical potential gradient into Fick's law combined with electro-osmotic drag. Here we additionally account for thermoosmosis by writing where λ is the hydration number of the ionomer, Vm its dry equivalent volume (Sec. 3.4), D λ the effective water diffusivity (Sec. 3.4), ξ the electro-osmotic drag coefficient (Sec. 3.5), and DT the thermo-osmotic transport coefficient (Sec. 3.6). The source term for dissolved water reads ω ∈ [0, 1] determines the mass fraction of water produced in dissolved rather than in liquid form at the ionomer-catalystcarbon triple phase boundary. Since the appropriate value for ω is an open problem [13,14], we assume ω = 1/2. The other half will appear as a source term for liquid water in Eq. 16. Assuming that (i) gas convection is negligible everywhere (∇P ≡ 0), (ii) no gas crossover through the membrane occurs, and that (iii) the supplied gases only consist of hydrogen on the anode side and an oxygen-nitrogen mix on the cathode side when dry, the total gas pressure P can be evaluated on these respective sides as The ideal gas law is assumed to hold, such that the partial pressures are given by pX = yX CRT , X = H2, O2, H2O, N2, where yX are the mole fractions of the gas species and C is the total interstitial gas concentration. Under these conditions and further assuming that thermal diffusion is negligible, an appropriate way to model the transport of gas species in PEMFC diffusion media is given by an in-series combination of the transport resistances of the Maxwell-Stefan model and Knudsen diffusion: where DX,Y denote the effective binary diffusivities and DK,X the effective Knudsen diffusivities in the porous layers (Sec. 3.8).
The inert nitrogen is not explicitly modeled on the cathode side because it follows from the requirement that X yX = 1. Species conservation is imposed with div jX = SX , X = H2, O2, H2O where the source terms account for reactant consumption, whereas accounts for the phase transitions from vapor to liquid or dissolved water and back. Liquid water transport is modeled with Darcy's law with the pore saturation s as the dependent variable, i.e, where Vw is the molar volume of liquid water, K abs the porous medium's intrinsic (absolute) hydraulic permeability, K rel the saturation-dependent relative permeability, ∂pc/∂s the layer's differential relationship between capillary pressure and saturation, and µ the dynamic viscosity of liquid water (all detailed in Sec. 3.9). To balance the evaporation/condensation sinks in Eq. 14 and to complete the electrochemical production of water halfway accounted for with Eq. 8, the liquid water source term must read

Phase change
Within the CLs, water is absorbed by the ionomer phase and desorbed again back to vapor in a relatively sluggish process (see Sec. 3.10). In order to account for this interfacial mass transfer resistance across the entire CL thickness, the sorption source terms in Eqs. 8 and 14 are modeled as where λeq is the equilibrium hydration number of the membrane (Sec. 3.7) and γ a,d = k a,d /L CL are the absorption and desorption rates with L CL the catalyst layer thickness (Sec. 3.13) and k a,d the ab-/desorption mass transfer coefficients (Sec. 3.10). Similarly, the phase change between vapor and liquid water is modeled as in which γe,c are evaporation/condensation rates (Sec. 3.11) and ysat = Psat/P denotes the saturation mole fraction of water vapor with Psat the saturation pressure (Sec. 3.7).

Contact resistance
Electrical contact resistance between individual components of a typical PEMFC is known to have a considerable stake in the overall fuel cell performance loss [15]. Here, we follow a recently demonstrated approach [16] to include both electrical (ECR) and thermal contact resistance (TCR) in the model. The continuity assumptions of the electric phase potential ϕe and temperature T at all affected layer boundaries are replaced by the constraints where [·] + ([·] − ) represents the value of the dependent variable [·] on the positive (negative) side of the interface as defined by the unit interface normal vector n, whereas Re (RT ) is the compression-dependent electrical (thermal) contact resistivity of the interface (see Sec. 3.13). ECR and TCR are modeled at the interior CL/GDL interfaces as well as at the exterior boundaries of the model, at the contact points between the GDLs and the bipolar plates (BPs), the latter of which are not an explicit part of the model. The PEM/CL interfaces are assumed to be perfectly thermally conductive and electrically insulating.

Boundary conditions
At the interfaces between the CLs and the membrane, zero-flux BCs are set for the electron flux as well as at the GDL/CL interfaces for the proton flux. The phase potentials ϕe and ϕp can freely be offset together, because only their derivatives and differences appear in the governing equations. We set ϕe = 0 at the anode current collector and impose the total cell voltage ϕe = U at the cathode current collector. Through electrical contact resistance (Eq. 19), this amounts to a coupling constraint between the electron flux and electron phase potential at each of the outer GDL surfaces. The temperatures of the anode and cathode bipolar plates, TA and TC, are assumed to be fixed. Water transport in dissolved form is bound to the ionomer phase. Therefore, zero-flux BCs are imposed on j λ at the GDL/CL interfaces. Since gas convection is neglected (∇P ≡ 0), the total gas pressure is given by the gas channel pressures PA,C. Dirichlet BCs are specified for the mole fractions of the gas species to match these conditions at the GDL/flow channel interfaces: yH 2 = 1 − yH 2 O at the outer AGDL surface and yO 2 = αO 2 (1 − yH 2 O) at the outer CGDL surface, where αO 2 denotes the mole fraction of oxygen in the dry supplied oxidation gas. At the same interfaces we additionally specify the relative gas humidity: yH 2 O = RHAPsat(TA)/PA at the outer AGDL surface and yH 2 O = RHCPsat(TC)/PC at the outer CGDL surface. No crossover of gas or liquid water through the membrane is considered, i.e., zero-flux BCs are set at both PEM/CL interfaces for all gas species and the liquid water saturation. At the CGDL/gas channel interface, the liquid water saturation is assumed to coincide with the immobile saturation sim (see Sec. 3.9).

Electrochemical reaction
The most common approach to model the reaction kinetics is the Butler-Volmer equation  [21] where S 0 A,C = j 0 A,C aA,C(1−s) is the product of exchange current density, reactive surface area density and a correction factor (1− s) for site blockage by liquid water. αA,C ( αA,C) are the forward (backward) half-reaction transfer coefficients. Measurement data for the hydrogen oxidation reaction suggests that they sum up to unity in the anode, at least for moderate current densities [17]. For the oxygen reduction reaction, on the other hand, there is no consensus in the literature on whether αC = αC = 1 holds, or whether there is a doubling of Tafel slope at intermediate voltages (αC = 1 to αC = 0.5) [18]. We use the former. The exchange current densities can be written as [17,18] where j 0,ref A,C are the exchange current densities at reference conditions (P ref = 1 atm and T ref = 80 • C), δA,C the kinetic reaction orders, EA,C the half-reaction activation energies, and finally, pA = pH 2 and pC = pO 2 are the reactant gas partial pressures in the two electrodes. The electrochemical parameters of the model are summarized in Tab. 1.

Electronic and ionic conductivities
Constant values are used for the effective electronic conductivities σ CL e = 390 S m −1 (for a catalyst layer with ionomer volume fraction i = 0.3) [22] and σ GDL e = 450 S m −1 (for a SGL 28 AA compressed by 1 MPa) [23].
The protonic conductivity of PFSA membranes and the ionomer phase in the catalyst layers has been the subject of more than 200 papers, with the majority focusing primarily on Nafion with equivalent weight (EW) 1100 [3]. For this reason, it is the first MEA material property that we pay particular attention to here, to examine the uncertainty in PEMFC performance prediction arising from scatter in experimental data on this material property. Being a strong function of the state of ionomer hydration, σp dominates the total ohmic losses in the cell in regions where the membrane is relatively dry. In Tab. 2, we summarize published parameterizations for Nafion 1100 EW, which can be (and have been) plugged into numerical models. Where an Arrhenius expression is used to account for temperature dependence, the reported activation energies are also given. Hsu et al. [24] originally carried a result over from percolation theory to express the protonic conductivity by a shifted power law, such that one can write where is the microstructure factor of the ionomer [40] ( i = τi = 1 in the PEM and i = 0.3, τi = 1.4 [41] in the CLs). The water RT [24] N117,  volume fraction in the hydrated ionomer is given by Later models by Thampan et al. [28] and Fimrite et al. [36] have extended the percolation-based conductivity by additionally expressing the prefactor σ0(T ) as a function of λ. Springer et al. [12], on the other hand, proposed a piecewise linear law in λ in their seminal modeling work -a simple correlation that is still widely used today. The phenomenological parameterizations used by Sone et al. [26], Yang et al. [34] and Maldonado et al. [38] use the water vapor activity (i.e., relative humidity) a to fit the observed ionic conductivity. To implement these into our model, the water activity in the bulk membrane is calculated by inverting the sorption isotherm λ(a) (see Sec. 3.7). Yang's relationship is excluded from our following analysis though, as it predicts far larger conductivities than all others and is considered to be an outlier for this reason. These proposed parameterizations of σp are juxtaposed in Fig. 1 to show that they vary considerably, although some spread may certainly partially be the result of different measurement temperatures. Despite characterizing the same class of PFSA membranes, the correlations by Hsu et al. [24], Morris & Sun [25], and Zhao et al. [39] predict far lower conductivity than the others. As noted by Weber et al. [3,32], the percolation model fits most measured data best. Their coefficients (f0 = 0.06, β = 1.

Thermal conductivity
The through-plane thermal conductivity k of the individual layers plays a subordinate role in MEA modeling, as will be shown in Part II. We therefore restrict the discussion on k to a single baseline parameterization for each layer type. In humidified Nafion membranes, it can be approximated by a linear interpolation between the thermal conductivities of water and the dry membrane [42]: : Review of protonic conductivities in vapor-equilibrated Nafion membranes. Below the largest λ at which σp(λ) = 0, the conductivity is set to vanish, which is omitted here for brevity. RT is short for room temperature. is a linear approximation to the thermal conductivity of dry Nafion, and is the internationally recommended correlation of the thermal conductivity of liquid water at 1 bar up to 110 • C with T = T /300 K [43]. A detailed review of available experimental data on the thermal conductivity of the remaining porous media can be found in [44]. Here, the Maxwell-Eucken equation [45] is used, where ks is the theoretical conductivity of the solid bulk material, k f the conductivity of the fluid filling the pore space, and p denotes the pore volume fraction. In order to obtain the bulk conductivity ks from effective conductivity measurements on real MEA materials, one can invert Eq. 28 using k f = 0 or k f ≈ 0.003 W m −1 K −1 , depending on whether the experiment was conducted in vacuum or air. For vacuum, this yields where k0 is the effectively measured thermal conductivity of the dry porous layer. Alhazmi et al. [46,47] have conducted the most comprehensive measurements of thermal conductivity of dry GDLs, k GDL 0 , as a function of temperature, clamping pressure and polytetrafluoroethylene content under vacuum conditions. We have fitted the following functional relationship to their data for the SGL 10 series, assuming that the correlations with temperature T and clamping pressure P cl are independent: where P ref = 1 bar. A power law was chosen for the pressure dependence because contact resistivities follow the same relationship (see Sec. 3.13). For the CLs a constant value is used for the dry thermal conductivity, k CL which lies approximately in the middle of the reported range of values in the literature. Humidity dependence is added through Eqs. 28 and 31 in the CLs just like in the GDLs. For the fluid conductivity k f , we assume that liquid water and the gas mixture form transport channels in through-plane direction along which heat is transported in parallel: This choice is motivated by measurements on humidified GDLs [49,50], which revealed an overall increase in effective thermal conductivity of as much as 50% at s = 0.25 and even more at higher saturations. The gas phase conductivity kg depends on the gas composition. We model it as a linear combination of the conductivities of the individual gas components, with the mole fraction as weights and species conductivities kX from [51,52]:

Water diffusivity in the ionomer
Water diffusion within and through the membrane has been a topic of extensive research over the past decades and has been thoroughly reviewed in the works of Kusoglu & Weber [2,3]. The numerous studies carried out to measure and parameterize the water diffusivity vary in the transport mode considered (driven by a pressure gradient, concentration gradient, temperature gradient, or chemical potential gradient), in the experimental technique (NMR, QENS, conductivity), in the types of diffusion coefficients being measured, in the studied regime (Fickean vs. non-Fickean, steady-state vs. transient), and even in the data separation and interpretation (e.g., correcting for measurement device resistance and interfacial resistance or not). As a result, the reported diffusion coefficients of dissolved water in Nafion are scattered over one to two orders of magnitude, which calls for a quantitative analysis of the modeling uncertainty associated with this material property. As shown in Fig. 2a,c, the intradiffusivity Di is a monotonically increasing function of λ [3], whereas most studies agree that the Fickean diffusivity DF exhibits a pronounced peak    [12] where none is given in the respective publication). All ordinates reach up to the self-diffusivity of water at 25 • C, 2.3×10 −5 cm 2 s −1 [74,75].   [77] 1250 EW [55] 1100 EW [58] N112 [78] N117 [61] Water self-diffusion [79] Parabolic fit (Eq. 35) Figure 3: Activation energy of the water diffusion coefficient in Nafion membranes as a function of membrane hydration.
around λ ≈ 2 − 4 (Fig. 2b,d). This local maximum in DF stems from the Darken factor, which relates the two diffusivities [12]: where λ(a) is the vapor uptake isotherm (see Sec. 3.7). Since λ(a) is relatively flat at λ ≈ 2 − 4, the Fickean diffusivity DF peaks there. However, the existence of a local maximum in DF is still a subject of controversy [63,73,76].
DF can either be measured directly or calculated from Di using Eq. 33. A compilation of published correlations between either of them and the level of hydration (via a, λ or fw) is given in Tab. 3 and Fig. 2c,d. Given DF, the effective diffusivity required for Eq. 7 can be calculated as with Mi from Eq. 23. We have recently identified the parametric expression given in the last row of Tab. 3 as the most plausible and convenient for numerical modeling of water transport within Nafion membranes [11]. It is used as the baseline parameterization also for the present study. For the temperature dependence, the Arrhenius equation is assumed to hold almost unanimously in experimental and modeling studies, even though a Speedy-Angell power law fits experimental data better [75]. The research group of Eisenberg [77,80] was the first to estimate the activation energy E d ≈ 20 kJ mol −1 for a Nafion membrane immersed in liquid water, which is close to the value for water self-diffusion (19.2 kJ mol −1 [79]). This value has subsequently been used by almost all researchers to fit their experimental data and to model DF(λ, T ), neglecting that it might depend on the mode and level of membrane hydration. Later measurements have brought up a variety of other values: 12-16 [60], 22 [73], 24 [61] and 28 kJ mol −1 [62]. Gong et al. [55], Kreuer et al. [58] and Guillermo et al. [78] investigated the dependency of E d on hydration. Their findings are juxtaposed in Fig. 3. A polynomial least-squares fit to the available data yields which is the activation energy implemented in the present model to adapt the rational polynomial in λ from Tab. 3 to arbitrary temperatures in the baseline parameterization. For a Nafion membrane of type NR211, which is considered for the baseline simulation, ρm = 1.97 g cm −3 [81] and mm = 1020 g mol −1 [82], such that the equivalent volume is given by   (22)) 5 [27] 2nd degree polynomial [84] bilinear [31] 2nd degree polynomial [33] 3rd degree polynomial [

Electro-osmosis
Electro-osmotic drag in Nafion has been another subject of controversial debate owing to the complexity of coupled ion/water transport and the difficulty in measuring it [3,4]. Springer et al. [12] proposed a linear approximation for the electro-osmotic drag coefficient based on a single data point for Nafion 117 immersed in liquid water: (ξ l , λ l ) = (2.5, 22). This relationship is still widely used in MEA modeling. Later measurements [87] have led researchers to conclude that vapor-equilibrated PFSA membranes are more appropriately characterized by ξ ≡ 1. Fig. 4 and Tab. 4 provide an overview over published measurements and parameterizations of the electro-osmotic drag coefficient, from which it becomes clear that a conclusive reliable correlation ξ(λ, T ) is still missing. Care must be taken when interpreting Fig. 4a, because some of these measurements were performed on liquid-equilibrated membranes with different degrees of drying to obtain mid-range water contents. These data might not be representative for vapor sorption. Agreement has not even been found on the general trend of the drag coefficient, with data suggesting an increasing, invariant, or even decreasing value of ξ with increasing membrane hydration.
A frequently used class of parameterizations are piecewise linear functions [31,102,103], which take the form when taking the theoretical upper bound ξ ≤ λ [83] into account.
Here, λ l is the water content for a liquid-equilibrated membrane, λv(1) the water content for a vapor-equilibrated membrane at unit vapor activity (see Sec. 3.7), and ξ l the drag coefficient for a liquid-equilibrated membrane. Eq. 37 is a simple attempt at taking Schroeder's paradox (λ l λv(1), see Sec. 3.7) into account in a global parameterization ξ(λ). Weber & Newman [32] suggested to use the Arrhenius equation with T ref = 30 • C and an activation energy of E ξ = 4 kJ mol −1 to model the temperature dependence of ξ l . In summary, Springer's linear law (Eq. 36) matches the widely scattered experimental data best. It is therefore is used here as the baseline parameterization, together with Eq. 38 for temperature dependence.

Thermo-osmosis
Thermo-osmotic transport of dissolved water is an entropic effect and occurs from cold to hot regions in hydrophilic membranes, i.e., in direction of the positive temperature gradient [104,105]. The thermo-osmotic transport coefficient DT is thus negative for typical PFSA membranes [106,107]. Kim & Mench [107] studied the temperature dependence of thermo-osmosis in liquid-equilibrated Nafion 112, Flemion SH50 and Gore-Select membranes and reported that the Arrhenius equation holds with an activation energy that is indistinguishable from that of concentration gradient-driven diffusion. This suggests that the transport mechanism for thermo-osmosis might be the same as for diffusion, albeit with different driving force. We therefore fitted their measured values of the thermo-osmotic transport coefficient in N112 using T ref = 80 • C and the activation energy E d from Eq. 35 and obtained Measurements of DT for vapor-equilibrated PFSA membranes are still missing in the literature. The magnitude of DT in various anion exchange membranes has been reported to increase with growing water content [108,109], and if the transport mechanism of thermo-osmosis is indeed similar to diffusion, it appears natural to assume that DT → 0 as λ → 0. Therefore, in order not to overestimate the effect of thermo-osmosis much in the model for vapor-equilibrated membranes, a linear approximation between zero and the reported coefficient for liquid-equilibrated Nafion is used here, analogous to Springer's linear interpolation of the electro-osmotic drag coefficient in Eq. 36:

Membrane hydration
In the parameterization of the equilibrium water content of the ionomer in Eq. 17, we account for simultaneous partial contact of the ionomer with liquid water and water vapor as well as for Schroeder's paradox by writing where λ l (λv) denotes the hydration number when the membrane is liquid-equilibrated (vapor-equilibrated). The vapor sorption isotherm λv(a) has been the subject of a vast number of experimental and theoretical studies. A chronological listing of proposed explicit and implicit functional relationships is given in Tab. 5. Historically, polynomial fits to experimental data (of degree 3 or higher in a) have been popular in PEMFC modeling, starting with Springer's and Hinatsu's measurements at 30 and 80 • C, respectively. Temperature dependence can be introduced by linearly interpolating between these two polynomials. In Fig. 5a, the different polynomials are plotted together with a selection of experimental data points from the open literature, showing that there is an uncertainty band of around 2 in width over the entire activity range, with a tendency to widen toward saturation. The water vapor activity is calculated as  [111,115,133], or multilayer adsorption theory, in its finite-layer Brunauer-Emmett-Teller (BET) form [28,116,134], in the infinite-layer limit [135], as a superposition of different sorption modes [126], or extended by elastic swelling [113,119]. Other swelling models Futerko & Hsing, 1999 [111] numerical solution of the implicit equation (1 + 0.2325a 2 (T /  use thermal, chemical and elastic equilibrium assumptions [121,123]. Many of these sorption models are based upon an additive decomposition of the total water content into a chemically bound (subscript b in Tab. 5) or clustered (c) part and a free (f) part.
Choi et al. [115] proposed a model in which they combine Flory-Huggins solution theory with elastic polymer swelling, using BET theory for the strongly bound water molecules. However, their description of the Flory-Huggins interaction parameter χ isn't explicit enough to allow us to reproduce their general sorption isotherm. Kreuer [123] proposed a thermodynamic sorption model with account for internal elastic pressure, which, after correction of the model equations [124], underestimates the membrane hydration due to the treatment of Nafion as an elastomer. These models are therefore excluded from our following uncertainty analysis.
In Fig. 5b the sorption models are plotted at 30 • C (where possible) to reveal their temperature-independent scatter, which is of an extent comparable to the polynomial correlations. Models that include a temperature dependence are also plotted in Fig. 5c at two different temperatures, highlighting that the effect of temperature is still far from understood. Four out of six models predict an increase in water uptake at high activity with rising temperature, whereas the other two predict the opposite. The experimental data these models have been validated on, as well as the proposed polynomial expressions, likewise disagree on the effect of temperature.
Thampan et al. [28] first recognized that finite-layer BET theory can be used to fit experimental vapor uptake data very well. Although Costamagna et al. [117] and Li et al. [122] later demonstrated that the Guggenheim-Anderson-de Boer (GAB) equation and a dual-mode (DM) model work just as well, the BET isotherm presents a good compromise between physical interpretation, high quality of fit, and suitability for model implementation (being an explicit relationship between a and λv, unlike the more complex implicit sorption models): where K denotes the ratio of the absorption equilibrium constant of the first layer to that of the subsequent layers, determining the shape of λv at low relative humidity. n is the number of adsorbed layers, governing the increase of water uptake at high relative humidity. The water loading at monolayer coverage λm can be estimated by [28] λm = a PEM p Aw Vm NA (45) where the area occupied by each adsorbed water molecule on the pore surface is approximately given by [136] NA is the Avogadro constant and a PEM p = 210 m 2 cm −3 the pore surface area per unit volume of the membrane [137]. Fitting Eq. 44 to published experiments on Nafion 117, NR211 and NR212 at temperatures between 20 • C and 30 • C [12,25,63,82,116,131,132] yields K = 92 and n = 12.8, which is close to the finding of Thampan et al. [28]. We use this BET isotherm for the baseline model parameterization. It also provides convenient oversaturation behavior, as depicted in Fig. 5d: As a → ∞, λv → λmn ≈ 23, which happens to coincide with the reported hydration number for liquid-equilibrated Nafion membranes. We therefore use λ l = λmn in Eqs. 36, 37, 40 and 41.

Gas diffusivity
An accurate way to estimate the effective binary diffusion coefficients DX,Y for Eq. 10 is given by Chapman-Enskog kinetic gas theory [52,138]. Assuming the ideal gas law, they can be calculated as where kB is the Boltzmann constant, mX denotes the molar mass of substance X, and is the microstructure factor of the pores, with compressiondependent porosity p and pore tortuosity τp as detailed in Sec. 3.13. For the saturation exponent φ, we use φ GDL = 3 [139] and φ CL = 1.5 [140]. A curve fit of the collision integral ΩX,Y is given by [141] ΩX where T * = kBT /εX,Y . For non-polar gases, the combining rules σX,Y = (σX + σY )/2 and εX,Y = √ εX εY can be used, where σX and εX are the Lennard-Jones collision diameters and potential depths, respectively, which are tabulated in the literature [138,142,143]. For inter-diffusion of a polar and a non-polar gas, modified theories with increased complexity exist, such as the one by Brokaw [144]. For simplicity it is assumed here that Chapman-Enskog theory also applies with sufficient accuracy for humid gases, i.e., that the dipole moment of water molecules can be neglected. Assuming cylindrical pores with effective average pore radius rp, the Knudsen diffusivities DK,X are given by [145] DK,X = 8rp 3 To account for pore narrowing by liquid water, we adopt the quadratic law rp = r p,dry s 2 w , which fits the data by Hutzenlaub et al. [146] well, where sw is the reduced wetting phase saturation as defined in Eq. 51. The following dry radii are used: r GDL p,dry = 15 µm for a SGL 24 BC [147] and r CL p,dry = 20 nm for a Nafion/carbon black CL [148].

Liquid water transport
The effective liquid water transport coefficient Ds in Eq. 15 depends on several material properties. In particular the functional relationship between saturation s and capillary pressure pc is a topic of extensive research with a large variety of fitting functions that have been proposed [149]. It strongly depends not only on the wettability and compression of the porous medium, but also on the exact kind of water transport process (primary injection, withdrawal, further injections). We found the overall impact of Ds on fuel cell performance to be mostly small, though, as will be shown in Part II of this series. Therefore, we focus on a single parameter set here rather than screening the literature for different parameterizations. A study on the effects of artificially altered capillary pressure-saturation relationships can be found in [150]. The van Genuchten law, which Gostick et al. [151] have found to apply to many common GDLs, is used here: with the following parameters for the secondary water injection curve of compressed SGL carbon paper: m = 0.6, l = 100, a breakthrough pressure of p b = 1.07 bar, and an immobile saturation of sim = 0.08 [151]. Mualem's model is used for the relative hydraulic permeability, reading [152] K rel = (1 − sw) 2 1 − s 1/m w 2m + 10 −6 (52) where the small offset serves to bypass numerical difficulties under dry conditions, i.e., to avoid that K rel → 0 as s → sim. We note, however, that the relative permeability in partially saturated GDLs is an active area of research, and alternate expressions such as S-shaped functions were recently proposed [153]. For the absolute permeability of the GDLs, the semiheuristic Carman-Kozeny equation for fibrous porous media is employed, as it proved to work well for carbon paper [154]:  Figure 6: Experimentally determined interfacial mass transfer coefficients for Nafion membranes as a function of water content. Individual symbols represent measurement data and lines show fitted functions. Data from [2,162] were converted from relative humidity to water content using Eq. 44.
with T = T /300 K. In the molar volume of liquid water Vw = mw/ρw, mw = 18.015 g mol −1 is the molar mass of water, and ρw its mass density given at standard atmospheric pressure by [128] ρw = ρc 1 + 1.9927 with T and ρc like in Eq. 43.

Vapor sorption
Numerous experiments have been carried out to measure the interfacial mass transfer coefficients ka and k d in Eq. 17, as listed in Tab. 6. They differ in the experimental procedure (absorption, desorption, permeation, liquid vs. vapor basins) and the driving force considered (concentration or hydration number gradient vs. activity gradient). Measurements that are based upon a water activity difference between the ionomer (ai) and vapor phases (a) are marked with an asterisk in Tab. 6. These coefficients can be approximately converted to ka and k d for use in Eq. 17 by dividing them by the local slope of the sorption isotherm, ∂λv/∂a [125], because for small deviations from equilibrium, Some of the more systematic studies that include the moisture dependence of the water transport resistance at the ionomer-gas interface are compared in Fig. 6. Ye & LeVan [72] reported an approximate power-law correlation with the partial pressure of water vapor. Since it is unclear how their expression translates to Eq. 17 though, we exclude this early result from the following analysis. Ge et al. [69] proposed to take ka and k d proportional to the ionomer's water volume fraction fw analogous to the intra-diffusion coefficient of water in bulk Nafion as proposed by Weber & Newman [32] (cf. Tab. 3). They also adopted the activation energy E ad = 20 kJ mol −1 from a diffusion measurement by Yeo & Eisenberg [80]. Later measurements of the activation energy of sorption have yielded larger values, E ad ≈ 29 kJ mol −1 [158,162,165]. He et al. [163] and Kusoglu et al. [2] proposed power laws in λ to express the mass transfer coefficients, based on data that possibly stem from the same measurement carried out at the Lawrence Berkeley National Laboratory, but assuming different residual hydration numbers λ0 (3.17 vs. 0 in Fig. 6). Converting their measured data points from activity-driven to hydration number-driven using Eq. 56 yields an approximately constant mass transfer coefficient k a,d ≈ 10 −3 cm s −1 independent of λ (data not shown).
Kongkanand [162] used a polynomial in water activity to fit the coefficients.
Although there is disagreement on whether the interfacial mass transfer grows with increasing λ [69,72,162] or not [2,163], the studies that differentiate between absorption and desorption agree that the former is substantially slower than the latter [69,125,162]. To date, the true dependence on moisture is essentially an open problem. Most researchers report the mass transfer coefficients in the range 10 −4 -10 −2 cm s −1 in the relevant temperature range of 50-90 • C. This is also the range numerically examined in early models by Okada [99,166]. With these considerations in mind, we adopt Ge's correlation for the baseline simulation.

Evaporation and condensation
Eq. 18 is the commonly employed out-of-equilibrium way to account for liquid-vapor phase change in macro-homogeneous PEMFC modeling with explicit representation of both phases. A comprehensive comparison of different expressions for the rates γe,c developed in modeling is given in Tab. 7. In the final column, they are evaluated for a typical state of fuel cell operation, demonstrating that the effective evaporation and condensation rates used in published two-phase models vary over as much as five orders of magnitude. Fig. 7 shows the timeline of when these rates came about in the literature, revealing no sign of convergence nor even a trend.
A feature shared by many models is that the condensation (evaporation) rate is assumed to be proportional to 1 − s (s) to account for the change of available phase boundary, as first proposed by Nguyen [168]. He et al. [170] proposed that condensation be proportional to the molar water vapor fraction yH 2 O, but more recently developed models have abandoned this assumption. Nam & Kaviany [174] suggested to use the Hertz-Knudsen equation from kinetic gas theory to model both directions of the phase change, assuming continuity of the temperature at the phase boundary and equal rates. Weber et al. [175] used a simpler expression, also based on equal rates for evaporation and condensation. These formulas explicitly include the liquid-gas interfacial area density a lg , which is unknown a priori, making even an estimation of the rates difficult. Nguyen [168] Song et al. [169] He et al. [170] Meng [171] Nguyen & He [172] Weber et al. [175] Birgersson et al. [176] Nam & Kaviany [174] Eikerling et al. [177] Wu et al. [178] Zenyuk et al. [179] Safi et al. [180] Lal et al. [181] [179,180] and Toray TGP-H-060 [181] were recast into the form of Eq. 18 and shown as heart symbols.
Wu et al. [178] estimate it to be a lg = Γsap(1 − s) for condensation and a lg = Γsaps for evaporation with an interfacial area accommodation coefficient Γs 0.2. ap denotes the average pore surface area density of the porous domain. A recent experimental study [179] has shown that a lg ∼ s and that there is no clear correlation between the evaporation rate per interfacial area and s, suggesting that the evaporation rate should indeed increase linearly with s.
It should be noted, however, that the general validity of the Hertz-Knudsen equation was recently questioned and that it was hypothesized that evaporation/condensation be governed by a balance of tiny pressure differences and momentum flux instead [182]. The Hertz-Knudsen rates used here are among the higher ones used in modeling works (cf. Fig. 7), quickly yielding saturated vapor where liquid water is present. In this quasi-equilibrium regime, the exact values of the rates become irrelevant [183,184], and any functional dependencies of the rate expressions on s, T etc. can essentially be dropped. This is in agreement with the reported observation that evaporation in partially saturated GDLs is fast enough to be transportlimited (diffusion-limited) rather than area-or kinetics-limited [179][180][181].

Latent heat
The molar latent heat of water condensation Hec used in Eq. 6 can be parameterized as This is a global least-squares fit to the tabulated data from [185]. The latent heat of water sorption can be written as where the mixing enthalpy Hmix, which is often assumed constant or even neglected in published models, is a function of membrane hydration [186,187] and temperature [188] for Nafion. Since a reliable parameterization for Hmix is absent from the literature, a new one is proposed here, based on the measurements on Nafion 115 by Wadsö & Jannasch [188]. Their data indicate that Hmix essentially decays exponentially in λ, with a temperature-dependent deviation in the very dry regime (λ < 3). We fitted the phenomenological approximation with temperature-dependent coefficients to their data in the range where T = T /300 K.

Compression and contact resistivities
To account for the effects of clamping pressure on the computational domain, the layer thicknesses are modeled as a function of pressure. Given the compressive strain εc, one can write where L0 denotes the thickness of the uncompressed layer. We model SGL 24 GDLs on both sides of the MEA, for which L GDL 0 = 190 µm [189]. The relationship between strain and applied clamping pressure P cl for these GDLs is [190] ε GDL c = −0.0083 P cl 1 MPa The two catalyst layers are also compressed under applied pressure, but reliable measurement data is rare. We fitted the following two-parameter function to the compressive strain of the CL reported by Burheim et al. [191]:  [196] 2.89 0.64 [197] Together with Eq. 61, this equation is used to express L CL as a function of applied pressure. For the baseline simulation, L CL 0 = 10 µm and P cl = 1 MPa are used. Nafion NR211 with constant thickness L PEM 0 = 25.4 µm [81] is chosen as the membrane, assuming for simplicity that swelling and compression cancel one another (ε PEM c = 0). Kumbur et al. [190] measured the average specific pore surface ap at three different compaction pressures from 0 to 1.4 MPa for a few GDLs from SGL Group with MPLs, finding a moderate increase with increasing pressure. Asymptotically, though, it is clear that at very large applied pressures, the pore surface area density must come down again. We therefore fitted a quadratic polynomial to their data for SGL 24 BC: where GDL p,0 = 0.75 [190] is the porosity of the uncompressed GDL. While Eq. 63 determines how much CLs are compressed, it is unclear how this changes their porosity. We assume that only the pore space of the CLs is compressed, which yields where CL p,0 is the porosity of the uncompressed CL, taken as 0.4 [192]. The relationship between trough-plane pore tortuosity and applied clamping pressure is commonly expressed indirectly, with tortuosity as a function of porosity. For the present model, a parabolic least-squares fit to experimental data for SGL 24 DA by Flückiger et al. [193] is used, whereas the tortuosity of the CLs is assumed constant (τ CL p = 1.5 [41,194]) due to the apparent absence of published measurement data on its compression dependence.
To complete the constitutive parameterization of our model, the electrical and thermal contact resistivities Re and RT in Eq. 19 remain to be specified for the different MEA interfaces. As recently highlighted [16], they usually follow power laws of the form where the coefficients R0 and ζ vary with the materials considered -a result that was also found theoretically for contacting fractal surfaces, for which the exponent is ζ ∈ [0.5, 1] depending on the degree of material plasticity and surface roughness [198]. Numerical simulations based on contact mechanics have confirmed this relationship [199,200] for ECR at the GDL/BP interface. Thermal and electrical contact resistances are implemented in the present MEA model using the values listed in Tab. 8. Finally, we note that through-plane TCR is significantly reduced by the presence of liquid water (e.g., [50]). Since the data in the literature is too scattered though to allow for a reliable parameterization w.r.t. s, the dependency of the contact resistivities on moisture is neglected here. This aspect requires additional experimental clarification.

Impact on fuel cell performance
In the model parameterization in Sec. 3, special attention was paid to six of the most critical and controversial transport parameters for Nafion-based MEAs, and their openly available constitutive relationships were reviewed in Tabs. 2 to 7. We now employ our macro-homogeneous steady-state MEA model to determine the degree of uncertainty in the fuel cell performance prediction associated with these six material parameterizations. Implemented in COMSOL Multiphysics, the model is numerically solved with the finite element method. 50 elements with quadratic Lagrangian shape functions are used per MEA layer (totaling in 250 finite elements for the entire MEA), and the damped Newton method is used with a relative error tolerance of 10 −5 to solve the coupled nonlinear differential equations simultaneously. Tab. 9 lists the operating conditions at which all simulations shown in this first part are carried out. They were chosen to approximately represent a differential section of a PEMFC operated in an automotive scenario. Fig. 8 shows the resulting polarization and power density curves predicted by the model when each of the six major parameters are substituted while all others are retained at the baseline as described in Sec. 3. As can be recognized from Fig. 8a, employing different expressions for the ionic conductivity σp of the Nafion membrane leads to enormous scatter. In order to quantify the discrepancy between the different model outputs, we selected the maximum reachable current density Imax and the peak power density Pmax as key figures. From the experimental scatter of σp alone, the total bandwidth of obtained values is 1.10 A cm −2 for Imax and 0.52 W cm −2 for Pmax, making the protonic conductivity of the ionomer a material property of very large uncertainty.

Protonic conductivity
The parameterizations by Hsu et al. [24], Morris & Sun [25] and Zhao et al. [39], which were shown to predict generally low conductivities in Fig. 1, yield current densities that are roughly a factor of two poorer than the others, over the entire range of cell voltages. Yet, if these three are dismissed as possible measurement outliers, the remaining performance variation is still appreciable (up to about 0.4 A cm −2 ). To reach large current densities, the ionomer needs to offer high conductivity at low water content, because electro-osmotic drag dries out the anode side of the membrane (cf. [11]). It is for this reason that the conductivity data by Maldonado et al. [38] yields high performance at intermediate voltages, followed by a sudden drop. The parameterizations which suggest higher σp at low λ values (in particular, those by Springer et al. [12], Costamagna [29] and Meier & Eigenberger [33]) outperform others at large current   densities. It must be stressed that even small differences in the conductivity expression that may appear marginal at first can have a great impact on fuel cell models and that detailed knowledge of σp(λ) at low λ is critical. Fig. 8b shows the influence of the water diffusion coefficient D λ . Generally, higher diffusivities result in stronger back diffusion, which helps keeping the anode humidified at large current densities, but the effect of nonlinear features in D λ can be subtle. Ye & LeVan's data [72] yields the poorest cell performance due to the very low diffusivity in the dry and wet regimes, despite being the largest in between (Fig. 2d). The correlation proposed by Lokkiluoto & Gasik [71], which globally predicts low diffusivity, also yields a low limiting current density. To maintain high performance at large currents, good water diffusivity under very dry conditions is crucial, similar to our conclusion for the ionic conductivity above. Fuller's [67], Kulikovsky's [31], Zhao's [62] and Caulk's [73] expressions yield the lowest membrane resistivites because their back diffusion coefficients remain the largest toward low λ values (cf. Fig. 2). Our own fit [11] to measurement data by Mittestadt & Staser [66] (with Arrhenius correction using Eq. 35) lies somewhat in the middle of the overall scatter of polarization curves. With 1.34 A cm −2 for Imax and 0.52 W cm −2 for Pmax, the uncertainty associated with D λ is even bigger than with σp.

Electro-osmosis
Next, we turn our attention to the electro-osmotic drag coefficient ξ. Generally, the fewer water molecules are dragged along with each traversing proton, the more evenly the membrane remains humidified across its thickness, which in turn reduces the ohmic resistivity of the membrane due to the monotonicity of σp in λ. It thus comes with no surprise that the large variation in measurement data on ξ translates to wide scatter in the resulting polarization curves. As Fig. 8c shows, the empirical constitutive expressions by Springer et al. [12], Dutta et al. [84] and Ge et al. [86] yield the best performance, since they predict the lowest ξ at medium to dry conditions (cf. Fig. 4b). It so happens, though, that this is also the regime where the reported experiments diverge the most (cf. Fig. 4a). When the parameterizations of Fuller [67,83], Eikerling [27], Meier & Eigenberger [33] or Lokkiluoto & Gasik [71] are employed, the fuel cell model stalls at much lower current densities, because the anode dries out more. In summary, the electro-osmotic drag coefficient in Nafion is a model parameter with large uncertainty.

Membrane hydration
A natural implication of the monotonic increase of membrane conductivity with increasing hydration is that higher sorption isotherms yield better performance prediction. This effect is apparent in Fig. 8d. Takata's [116] and Costamagna's [117] vapor sorption models are among those which yield the smallest λv values in the low and high activity regimes, respectively (cf. Fig. 5). They therefore result in the steepest decline of the polarization curve in the ohmic region. Morin's data [127] and Kusoglu's model, on the other hand, both suggest good membrane hydration over the entire activity range, resulting in a polarization curve that extends toward higher current densities. With a total spread of 0.60 A cm −2 for the tested parameterizations, the water vapor uptake of the membrane is a significant source of modeling uncertainty. Given that it depends also on the membrane's hygro-thermal history [38,201] (which has been ignored in our present analysis), this highlights that detailed experimental characterization of the membrane is required to make PEMFC models predictive.

Vapor sorption kinetics
The next effect in the line is the sorption of water vapor at the ionomer-gas interface. With the four hydration-dependent expressions for the mass transfer coefficients plotted in Fig. 6, we estimate the modeling uncertainty originating from measurement data on k a,d , setting ka = k d for the two parameterizations that do not distinguish between absorption and desorption [2,163]. Fig. 8e shows the resulting polarization curves, which start to separate no earlier than at intermediate cell voltages. He's expression [163] yields the most limited cell performance due to the fast-dropping mass transfer coefficients toward low hydration numbers. The slight kink near 0.6 V is due to the kink in k a,d at λ0 = 3.17, which is used here with the intention of adopting He's expression without modification, in spite of the absence of residual hydration in the employed sorption isotherm (Eq. 44). Konkanand's data [162] yields the largest current densities, for their mass transfer coefficients remain positive (still allowing for moderate vapor absorption) when the electro-osmotic drag dries out the anode side of the membrane. Once again, we find that the constitutive material behavior at very low water content has a major impact in the predicted cell performance and therefore needs to be known with high accuracy -higher than currently available in the literature.

Evaporation and condensation
Finally, the effect of adopting different expressions for the evaporation and condensation rates is shown in Fig. 8f. Wu's Hertz-Knudsen expression [178] is used as the baseline parameterization with lowered evaporation coefficient (see Sec. 3.11), which results in higher liquid water saturation, and consequently, lower limiting current density. With a difference of just 31 mA cm −2 , the effect is relatively small, though. Expressions which give fast evaporation and slow condensation, such as [168,170,176] (cf. Fig. 7), reduce the liquid water saturation further, in particular in the CCL, allowing for better oxygen access to the catalyst sites. In total, though, the scatter induced by the different phase change rates is much smaller than for the five ionomer properties discussed above, because in the examined parameter range, evaporation is fast enough to yield a relative gas humidity of ≈ 100% across almost the entire MEA where liquid water is present.

Result summary
The results of our uncertainty analysis are summarized in Tab. 10, where in addition to the total spread, also the standard deviations for the limiting current density and peak power density are given. This corroborates that the same ranking list is obtained when individual parameterizations, which may be deemed measurements outliers, are given less weight or omitted from the analysis. Evidently, getting the water transport inside the membrane and across its interfaces correct is the key to predictive macro-homogeneous MEA modeling for PEMFCs. Most uncertainty stems from the ionic conductivity as well as from the two major water transport effects in the ionomer: back diffusion and electro-osmosis.

Conclusion
With the present work, we have found a quantitative answer to the questions raised in the introductory section: When characterizing a MEA with the intention of extracting material properties for PEMFC modeling, special attention should be paid to the membrane properties, in particular the Fickean diffusivity of dissolved water, the protonic conductivity and the electro-osmotic drag coefficient. At least under the present modeling assumptions and operating conditions, these are the top three traits of Nafion which cause the most uncertainty in the predicted fuel cell performance, based on the available data in the open literature. For high accuracy at large current densities, their functional dependency on the state of membrane hydration is critical, in particular in the dry regime (low λ). This finding is relevant not only for fuel cell experimentalists, but also for modelers who are faced with the task of selecting appropriate MEA material parameterizations for their calculations.
With the present scatter in the available data records on Nafion properties, predictive performance prediction is very difficult. Measurement errors should be quantified more routinely in experimental studies than they have been in the past. Further research is required to improve the understanding of the transport processes of hydrogen ions and water molecules across the electrolyte membrane, not only under varying water content, but also to clarify the temperature dependence of water uptake and electro-osmotic drag.
In Part I, we have studied only local parameter sensitivity, i.e., at fixed operating conditions and with all material parameterizations but one fixed. Given the highly nonlinear nature of PEMFCs, the picture might change in different scenarios. So far, we have focused on uncertainty in the extrema of polarization characteristics, which is induced by scatter in the proposed parameterizations for six of the most critical MEA properties. In Part II, we perform a global parameter sensitivity analysis to gain insight into the general model response to changes in the MEA parameters, not only limited to cell performance, but also regarding the predicted heat and water balance.