A Mixed Wettability Pore Size Distribution Based Mathematical Model for Analyzing Two-Phase Flow in Porous Electrodes Mathematical Model Improving the performance of polymer electrolyte

A multi-dimensional, non-isothermal, two-phase membrane electrode assembly (MEA) numerical model is developed where the micro-structure of the porous layers is characterized by a mixed wettability pore size distribution (PSD). The PSD model is used to predict local water saturation based on gas and liquid pressures, and can be used to study the effect of varying pore size and wettability. The MEA model accounts for gas transport via molecular and Knudsen diffusion, liquid water transport, sorbed water transport by back-diffusion, electro- and thermo-osmosis, and heat generation and transport. Multi-step kinetic models are used to predict anode and cathode electrochemical reactions. Local transport losses are accounted for using a local transport resistance. The PSD model is used to predict capillary pressure vs. saturation and saturation vs. relative liquid permeability curves based on PSDs from several GDLs obtained using mercury intrusion porosimetry. The PSD-based MEA model electrochemical performance predictions are also compared to experimental data from the literature. Results show that the numerical model is able to capture the performance changes associated with varying temperate and the introduction of a micro-porous layer.

Improving the performance of polymer electrolyte fuel cells (PE-FCs) at high current density is critical for reducing stack size, cost and weight in automobile applications. Water accumulation in the electrode however, limits fuel cell performance at high current densities. [1][2][3] In order to mitigate the excessive water buildup in the cell, which hinders mass transport, gas diffusion layer (GDL) and catalyst layer (CL) microstructures and wettabilities must be advanced to achieve sufficient water retention in the electrolyte, reject excess water in vapor form and alleviate complete flooding of the electrode. For example, by modifying the weight percentage (wt%) of the hydrophobic content in the GDL, Wang et al. 4 reported a 10%-20% increase in performance. To understand fuel cell behavior in the two-phase region and optimize the design of GDL and CL, a two-phase model which includes microstructural information of the porous layers is needed.
Two approaches have commonly been used to study liquid water transport in fuel cells. The first approach involves the solution of a saturation equation obtained by replacing liquid pressure by saturation in Darcy's law. [5][6][7] In this approach the Levertt J-function is commonly used to relate saturation with capillary pressure. The Leverett J-function is however, an empirically measured curve for a given material and it is obtained from porous media structure analysis. Thus, it is difficult to optimize the structure of the porous layers using Leverett J-functions. The second approach, known as the multi-phase approach, involves solving mass and momentum equations for gas and liquid phases independently and then relating the capillary pressure to saturation using a closure equation based on micro-structural information such as a pore size distribution (PSD). [8][9][10][11][12] The use of gas and liquid pressure naturally results in continuous solution variables at the interface between layers with different wettabilities which is not usually the case with the first approach based on saturation. [5][6][7] Models relating gas and liquid pressure to saturation by means of a PSD have been developed in recent years. [10][11][12][13] In order to develop a multi-phase model, a set of closure equations that use the PSD as a metric to characterize the material and estimate transport properties have been developed. The first structure-based models that proposed the use of a PSD assumed either a hydrophilic or hydrophobic pore network for the whole material. 10,11 More recently, researchers have introduced a variety of mixed wettability models based on experimental observations showing that a network of hydrophilic and hydrophobic pores z E-mail: secanell@ualberta.ca might co-exist in fuel cell materials such as GDLs and CLs, e.g., Refs. 14, 15 Weber et al., 11,16 Mateo et al. 13 and Mulone et al. 12 have treated hydrophilic (HI) and hydrophobic (HO) pore-networks separately either by introducing a continuous wettability distribution, i.e. Ref. 16, or by studying two independent networks, e.g. Refs. 11,12 and 13. In all previous work, due to the difficulty of implementation of a PSD model in multi-dimensional computational fluid dynamic solvers, the PSD model was either not integrated in a complete membrane electrode assembly (MEA) model 10,12,13 or was integrated only in a one-dimensional model. 11,16 A two-dimensional model is however required in order to study the water distribution in channel and land regions. Water vapor fluxes have also seldom been studied in two-phase models in the literature. Finally, previous models were only shown to operate under cold and wet conditions and its suitability under dry conditions was not demonstrated.
Evaporation and condensation of water in different regions of the cell leads to phase change induced (PCI) flow. PCI flow has been studied in the literature either by experimental in situ measurements 17,18 or using 1D non-isothermal, two-phase models to simulate water evaporation/condensation quantitatively. 19 The temperature-gradientinduced water flux measurements conducted by Mench et al. 17,18 through a non-operating and operating fuel cell show the importance of the PCI flow at various temperatures. Weber et al. 19 coupled thermal and water management together with a non-isothermal 1D model to investigate the transport of water due to the temperature gradient. However, a 1D model cannot fully capture channel/rib effect on water movement due to non-uniform temperature distribution across the MEA as observed by Bhaiya et al. 20 A detailed multi-dimensional analysis of the water fluxes in liquid and vapor form and the areas of evaporation is provided in this article with the aim at helping to understand the importance of PCI flow and to study under which conditions water is evacuated in liquid or vapor form.
In this article, a two-dimensional, non-isothermal, two-phase flow model that takes into account the microstructural information of the porous layers in the MEA is developed. The model can be used to study the effect of the PSD and the mixed wettability of different layers on fuel cell performance. It is implemented in the open-source program OpenFCST. 21 The structure of the article is as follows. First, the MEA mathematical model, the numerical PSD closure equations, solution strategy, post-processing analysis and input parameters are discussed. Then, the methods used for measuring the PSD are outlined. In the last section, the PSD model liquid transport property predictions as well as the two-phase model electrochemical performance predictions are validated by comparing to literature data.

Mathematical Model
Assumptions.-The mathematical model used to study two-phase flow transport is subject to the following assumptions: (1) The fuel cell is at steady-state.
(2) Gas convection in the porous media is neglected as the gas pressure is assumed constant. (3) The hydrophilic and hydrophobic pore-networks are homogeneous in the diffusion medium in each representative element volume (REV). (4) Liquid water movement is governed by capillary driven flow. (5) The liquid water is considered an incompressible and Newtonian fluid in the laminar regime in the porous media. (6) The phase change in the pore-network is driven by the gradient between partial pressure of water vapor and saturation water vapor pressure. (7) Membrane is impermeable to liquid water. Assumption 3 is valid in the CL but is only an approximation in the GDL. The latter is justified on the basis of the averaging performed along the channel to create a two dimensional model. Assumption 4 is based on the dimensional analysis in Ref. 24 for liquid water in the fuel cell. The capillary number Ca is found to be of the order of 10 −7 , Bond number is 10 −3 , and Weber number is 10 −10 . Hence, viscous drag, gravity effects and inertial effects are negligible compared to surface tension in this case. Assumption 6 is commonly used in the literature 10,11,13 which employs a single equation to account for both condensation and evaporation.
Governing equations.-The non-isothermal MEA model proposed by Bhaiya et al. 20 and implemented in OpenFCST was extended to include two-phase flow. Darcy's law is used to predict liquid pressure such that the overall system of equations is: where c tot is the total concentration of the gas mixture, x i is the molar fraction of species i, D eff i, j is the effective diffusivity of species i in j (based on molecular and Knudsen diffusivity), φ m and φ s are electrolyte and electronic potentials respectively, σ eff m and σ eff s are the corresponding effective conductivities, λ is the ionomer water content, D eff λ is the sorbed water diffusion coefficient, and J i andH i are the molar diffusive flux and the specific enthalpy of species i. The last equation describes the transport of liquid water where k lr is the effective liquid permeability, μ l is water dynamic viscosity and ρ l is the density of water. The capillary pressure is defined as: where p l is liquid water pressure and p g is the pressure of air and it is assumed to be constant. The source or sink terms corresponding to each equation are given in Table I. The terms S heat andẆ electrical are detailed in the previous non-isothermal model. 20 j is the volumetric current density computed using the recently proposed ionomer coated catalyst particle (ICCP) model. 25,26 Evaporation or condensation is considered to take place at the liquid-vapor interface such that the term S H2O is given by: , otherwise [3] where k e and k c are the volumetric evaporation and condensation rate constants in mol · cm −2 s −1 , a lv is the liquid-gas interfacial surface area per unit volume in cm 2 /cm 3 , and p v is the vapor pressure. The effective saturated vapor pressure in a capillary, p sat K( pc ,T ) , is determined by considering the Kelvin effect and the Young-Laplace equation as follows: where p sat (T ) is the uncorrected saturated vapor pressure of water, M H 2O is the molar mass of water, ρ g is the density of water vapor, R g is the universal gas constant, and T is the temperature. The additional source or sink term added to the thermal equation corresponds to the latent heat of vaporization or condensation expressed in the form of: [5] where H lv is the exchange molar latent heat in vaporization or condensation of water in kJ/mol. In the thermal equation, an additional term is added to account for the enthalpy transport of liquid water: where J l is the liquid water flux and H l is the molar enthalpy of liquid water.
Closure equations.-Closure equations are used to relate liquid and gas pressures and layer microstructure to saturation and other effective properties. Using the commonly used bundle of capillaries F532 Journal of The Electrochemical Society, 164 (6) F530-F539 (2017) idealization of the porous media, 10,11,13 pressures and saturation are related via a pore size distribution (PSD).
Pore size distribution.-The PSD is divided into two individual lognormal distributions over the pore sizes representing the hydrophilic and hydrophobic pore networks respectively. The initial form of the mixed wettability PSD model is: where V T is the total pore volume and dV (r )/dr is the differential accumulated volume with respect to a change in pore radius dr. X (r ) H I is the normalized volume of hydrophilic pore per total pore volume. The expanded PSD model takes the form: where E HI,k and E HO,k are: F HI and F HO are the volume fraction for hydrophilic and hydrophobic pores respectively, f i,k is the partial contribution of sub-distribution functions from mode k in either hydrophilic or hydrophobic network (referred here with index i), r i,k is the characteristic pore size of subdistribution from mode k, and s i,k is the standard deviation.
Saturation.-Saturation is estimated by integrating the contributions from all pores in the porous media as follows: where r c is the critical radius estimated using the capillary pressure and Washburn equation. 27 The limit of integration takes into account that, in a hydrophilic media, small pores would be invaded by water first, followed by the larger pores. The opposite occurs in a hydrophobic material. Analytically solving the integral above, an explicit equation for saturation is obtained 13 : Permeability.-Absolute permeability.-The absolute permeability is calculated using Hagen-Poiseuille equation and integrating the contributions from all pores in the layer: where ε o is the porosity and λ PSD represents the interconnectivity of the pores. The term in front of the integral represents the probability of finding two interconnected pores. Combining the probability with the previous Equation 11 and solving the integrals, the absolute permeability in a fully saturated porous medium is: Relative liquid permeability.-The relative liquid permeability is determined as the ratio between the effective liquid permeability and the absolute permeability as follows: [15] where all terms except k sat accounts for the contributions of a partially saturated layer from both hydrophilic and hydrophobic pores.
The integration forms to get hydrophilic and hydrophobic contributions can be found in the Appendix.
Liquid-gas interfacial surface area.-To account for the possibility of having a common area for liquid and gas between filled and empty capillaries, a probability density function is defined: where a(r ) c is the liquid invaded cross-section area per unit volume of cylindrical pores and a max is the total cross-section area per unit volume which is estimated as: The probability function forces the liquid-gas interfacial surface area to be zero either when the layer is fully dry or saturated, whereas when the layer is partially saturated, there exists a maximum value for the interfacial surface area. The cross sectional area per unit volume is determined by taking the ratio between the cross-section area (πr 2 ) of the pore and its volume (πr 2 L): [18] where V T is the total pore volume and L is the length of the cylindrical pores which is estimated as 4r as proposed in Ref. 10 due to the lack of experimental data. The interfacial surface area expressions for hydrophilic and hydrophobic pores can be found in the appendix A. The overall expression for the interfacial surface area per unit volume is: [19] Diffusion.-Average Knudsen Radius.-The average capillary radius is estimated by taking the ratio between pore volume and the lateral area surface, and multiplying by two to get the radius:  [20] where the lateral area per unit volume is given by:  [21] The radius calculated for Knudsen diffusion accounts for the collision between gas molecules and pore walls. Therefore, the limits of the integral over the pores are switched to estimate the volume of gas pores. The detailed integral expressions are in the appendix A. Diffusivity.-Inside the pore space of the CL, the length scale of the capillary radius is comparable to the mean free path of the diffusion particles and Knudsen diffusion is important. In the GDL, where the pores are larger, molecular diffusion is dominant. The Knudsen diffusion coefficient of species i is estimated as: where M i is the molar mass of species i, T is the temperature and R g is the gas constant. The effects of molecular and Knudsen diffusivity are combined using the Bosanquet equation 28 : where D ij is the binary molecular diffusion coefficient. For both GDL and CL, the effective diffusivity is a function of saturaiton. In the GDL, the partially-saturated carbon fiber diffusion layer model proposed by Garcia-Salaberri et al. 29 is used for both in plane and through plane directions. The effective gas diffusivity is: [24] where is the porosity, τ is the tortuosity, s is the local saturation and n is the exponent which characterizes the impact of liquid water saturation on diffusive gas transport, i.e., n is 3.5 in the through-plane direction and n is 2.5 in the in-plane direction. The effective diffusivity of the CL is estimated using percolation theory 30 : where cl V is the volume fraction of void space in CL, γ is the exponent which accounts for the loss of mass transport due to liquid water, th and μ are constants that depend on the orientation of the components in CL and the function θ( V (1 − s) − th ) is the Heaviside unit step function. In this article, μ is 2.0 based on Ref. 31, γ is assumed to be the same as μ, and th is 0.25 as it provides a good agreement with the experimental data reported by Yu et al. 32 as shown in Figure 2. Figure  2 shows the predicted value using the expression in the article for the CL under study (r Kn = 100 nm) as well as smaller Knudsen radii. The smaller Knudsen radii in Figure 2 indicates that when the average Knudsen radius is reduced under saturated conditions, the effective oxygen diffusivity might be strongly influenced by Knudsen effect.
Boundary conditions.-The two-dimensional computational domain is illustrated in Figure 1. The boundary conditions used for oxygen, water vapor and temperature are given in Ref. 20. The selection of appropriate boundary conditions at the GDL-channel interface is critical to any two-phase flow model. Zero flux 7,9,33 and fixed saturation values are usually imposed as boundary conditions. 34 These boundary conditions are either forcing all water to be vaporized in the MEA or imposing the existence of liquid water at the GDL-channel boundary. Their applicability is therefore limited and a more general boundary condition should be developed. An example of a more appropriate boundary condition is that recently proposed by Zenyuk et al. 2 where a step function is used. After reaching the breakthrough pressure, the liquid pressure however, is considered to be constant in this method. A dynamic boundary condition which relates liquid pressure to liquid water flux might be a more realistic choice and is implemented in this article. If the capillary pressure is below a given threshold value, i.e., a breakthrough pressure, a zero liquid water flux is imposed. Once the capillary pressure reaches a given breakthrough pressure, 2000 Pa in this article based on Ref. 35, a flux proportional to the liquid pressure is applied, i.e., ρ l u l · n = − ρ l k rl μ l ∇ p l · n = k p l − p l,channel p 0 g( p l ) [26] where g( p l ) = tanh(( p l − p l,channel )/ p 0 ) + 1 2 θ( p l − p BT ) , [27] k is an unknown proportionally constant that controls the flux of water as a function of the liquid pressure, p BT is the liquid breakthrough pressure, and θ( p l − p BT ) is a step function, i.e., it is set to be zero until p l > p BT is satisfied in the Newton solver loop and not modified further in order to maintain numerical stability. Its validity is confirmed during post-processing by making sure that the liquid water flux remains positive. By performing a parametric study on k, a value of 10 −6 kg/(m 2 · s) is selected. In order to prevent negative liquid water flux entering from the boundary once the step function is set to be one, a tangent function tanh(( p l − p l,channel )/ p 0 ) is used. p l,channel is the liquid pressure at the channel-GDL interface and it is set to be atmospheric pressure considering the droplet volume is large enough so that Laplace pressure is negligible, and p 0 is a non-dimensional factor which is set to be 1 Pa in this case. No liquid flux boundary conditions are applied at the PEM/CL and GDL/land interfaces and at symmetric boundaries, i.e., upper and lower faces in Figure 1. It is assumed in the model that the dominant water transport mechanisms through the membrane are back diffusion, electro-osmotic drag and thermo-osmosis. Liquid water permeation through the membrane is not considered even though, under fully humidified conditions, membrane liquid water permeation has been observed. 36,37 Future work should aim at improving the membrane model to include the effect of liquid pressure.
Solution strategy.-The governing equations were implemented in the open source fuel cell simulation software OpenFCST. 21 The non-linear system of equations is linearized using Newton's method. The linearized second order partial differential equations are then discretized using the Galerkin weighted residuals method. Second order Lagrange elements provided by the deal.II libraries 38 are used to approximate all solution variables. The system matrix and the force vector are evaluated using Gaussian quadrature rule. The resulting non-symmetric matrix is solved using the multi-frontal parallel distributed solver (MUMPS). 39 The initial mesh contains 12,775 degrees of freedom. The mesh is adaptively refined once using the adaptive error estimator by Kelly et al. 40 in deal.II by dividing the 30% of the cells with the largest errors into four cell in order to minimize local errors. The convergence criteria for the nonlinear solver is that the normalized L2-norm of the residual is smaller than 10 −8 , i.e., The solutions are grid independent as the difference on integrated quantities, e.g., total current density and water flux at boundary, at each cell voltage between the current mesh and a systematically refined mesh is within the order of 10 −6 . The overall running time to achieve a complete polarization curve with 20,832 degrees of freedom with 4 processors in parallel is between 8.37 and 651.83 minutes, depending on convergence rate. Convergence is challenging in the regions where the cell transitions from dry (only water vapor leaving the cell) to wet (liquid water leaving the cell mode).
Post processing.-Post-processing routines were implemented in OpenFCST 21 to evaluate velocity vectors and volumetric and surface integrals based on the solution over the computational domain. Water fluxes are computed at post-processing by solving the following boundary integrals: The computed water fluxes are also used to confirm that the conservation equation is satisfied globally. The evaporated or condensed water in the cathode electrodes, current density and heat generated are computed by integrating the volumetric source terms in Table I over the cathode CL, e.g., The cell resistance which is usually measured either by current interrupt or high frequency impedance analysis can also be estimated in the post-process routine. 26 Using the calculated heat dissipated in the GDL, MPL and CL due to electronic transport and heat dissipated in the PEM due to protonic transport, the overall cell resistance is estimated as 26 : [31] whereq disp is the total heat dissipation, i is the current density and A active is the experimental active area. The model predicted cell resistance is likely lower than the actual cell resistance since contact resistances are not accounted for.
Input parameters.-The MEA parameters, such as computational geometry, kinetic parameters and transport properties, are set such that it reflects the PEFC configuration. Channel and current collector width are 1.0 mm. Table II shows PSD GDL and CL baseline parameters in the simulations. The porosity and PSD of Toray 090 GDL (20%PT F E) and SGL 34BA are measured in-house using mercury intrusion porosimetry (MIP). The experimental data for Toray 090 is in agreement with a previous publication. 41 The volume fractions of hydrophilic pores in the GDLs are estimated by manually fitting the capillary-saturation curves in the literature. 35 The CL PSD is measured using MIP and the thickness is measured using scanning electron microscopy (SEM). The CL porosity is estimated by subtracting electrolyte and solid phase volume fractions from one. The external CL contact angle has been measured in references. 42,43,44 However, a direct measurement of the CL internal contact angle has not been reported. In this study, the hydrophobic contact angle used in the CL is 93 • and it is an average of the values (92 • and 94 • ) reported by Yu et al. 43 The contact angle assigned to the hydrophilic pore-network is 84 • by considering the contact angle of Ketjen Black carbons reported by Soboleva et al. 45 All electrochemical parameters for CLs are obtained from the experiments discussed in the MEA fabrication and testing section in a follow up article. 46 The thermal related parameters, i.e., GDLs, CLs and PEM thermal conductivity as well as all PEM transport properties used in the model are from the previous OpenFCST non-isothermal MEA model validation study reported by Bhaiya et al. 20 The double trap and dual path kinetic parameters used in the model are given in Moore et al. 22 and Wang et al., 47 respectively. The micro-scale ICCP model parameters used in the simulation were reported by Secanell et al., 26 in particular the local resistance rate constants are 0.1 and 0.001 cm/s in anode and cathode, respectively for 80 • C and 50%R H, and 60 • C. For 80 • C and 90%R H the values are 0.1 and 0.005 cm/s. In order to obtain the value for local resistance rate constant at 80 • C and 50%R H and 80 • C and 90%R H, simulations for CLs with varying Pt loading were performed using multiple local resistance rate values until the slope of the total oxygen resistance vs. the roughness factor (i.e., the product of the Pt loading and electrochemical active area) for the model and the experiments in Reference 48 were similar. The values of the local resistance rate that provided the best fit for the two relative humidities are then used in the study.
All input parameters in the model are the same for all simulations except evaporation and condensation rates which are considered to be temperature dependent, and the local resistance constant used in the ICCP model 26 at 80 • C and 90%R H. At 80 • C, 50%R H and 90%R H cases, an evaporation rate of 2 · 10 −3 mol/(cm 2 · s) is used as it allows most water to be evaporated, as observed experimentally. At 60 • C and 40 • C, a value of 2 · 10 −4 mol/(cm 2 · s) is used in order to take into account experimental evidence reported by Zenyuk et al. 49 showing a decrease in evaporation rate with temperature. At 80 • C, a condensation rate of 2 · 10 −3 mol/(cm 2 · s) is used as smaller values result in an RH greater than one. At 60 • C and 40 • C, the condensation rate constant used is 2 · 10 −2 mol/(cm 2 · s).

Experimental
Pore size distribution estimation.-Pore size distribution is determined using mercury intrusion porosimetry (PoreMaster 33, Quantachrome). For GDLs, two 5 cm 2 samples were cut into strips and inserted in the penetrometer. For CLs, three double sided 5 cm 2 MEA samples were used. The penetrometer was then evacuated up to an absolute pressure of 0.004 psi, and then, further evacuated for 30 minutes. The cell was then filled with mercury and pressurized to a maximum pressure of 33000 psi. The intruded volume of mercury and pressure were recorded. The Young-Laplace equation is used to estimate the PSD from the cumulative pore volume curve. The machine pore diameter range is between 950 and 0.0064 μm.

Results and Discussion
Validation of PSD-based closure relationships.-Gas diffusion layer.- Figure 3 shows the experimentally measured and numerically fitted PSDs for SGL 34BA and Toray 090 (20%PT F E). The method of least squares and Equation 5 are used to fit the experimental data in a LibreOffice spreadsheet. It can be seen from Figure 3 that the numerical fittings agree with the experimental PSDs in both cases. The pores in SGL 34BA (see Figure 3 (a)) are centered at 14.2 μm and 34.0 μm, respectively. It has a wider pore size range as compared to Toray 090 which is characterized by two modes centered at 18.0 μm and 23.0 μm (see Figure 3 (b)). The percentage of hydrophilic pores assigned to each GDL is obtained by fitting the residual liquid saturation. Although the residual saturation cannot be interpreted as evidence of hydrophilic pores, since it might also be due to pore geometry and connectivity, it is used in the model to represent the hydrophilic pores because it accounts for the liquid water that cannot be easily removed from the porous material. The smaller PSD mode is assigned to the hydrophilic phase as shown in Figure 3.
The capillary pressure vs. saturation relationships obtained from the PSD model are compared against the experimental data of several GDLs in Figure 4. First, to validate the PSD implementation, the experimentally measured capillary pressure-saturation relationship for SGL 10AA is reproduced using the parameters reported by Weber et al. 16 With the correct implementation, the capillary pressuresaturation relationship for SGL 34BA is then compared to the experimental data of SGL 10AA and SGL 10BA reported by Gostick et al. 35 The numerically predicted capillary pressure-saturation relationship for Toray 090 is also compared to the experimental data reported by Gostick et al. 35 and Fairweather et al. 50 in Figure 4b. The predicted values for SGL 34BA and Toray 090 are within the values reported for the same type of GDLs in Ref. 35,50. The model seems to provide reasonable results. It is known that the capillary pressure vs. saturation relationship shows hysteresis. 35,50 Only the imbibition branch is used in this article because it is assumed that water enters the layer following an imbibition process but is usually removed by evaporation instead of drainage in an operating fuel cell.
The relative liquid permeability predicted by the Toray 090 and SGL 34BA PSD models are shown in Figure 5. The values are between the experimental measurements reported by Luo et al. 51 and Koido et al. 34 The absolute permeabilities of the two GDLs have been matched with literature reported values 41,52 by setting the interconnectivity factor λ P SD in the PSD model to be 1.26 and 2.00 for SGL 34BA and Toray 090 respectively.
The bundle of capillaries idealization appears to provide reasonable agreements with the experimental data thus far. The connectivity between the pores is accounted for only when estimating the absolute permeability using λ P SD in Equation 11. The λ P SD reflects interconnection between randomly rejoined capillaries and is determined by fitting the experimentally measured absolute permeability. In estimating the rest of the parameters, i.e., saturation, liquid-gas interfacial surface area and average Knudsen radius, information about connectivity is lost in the PSD model as compared to a pore network 53 or micro-scale simulation. 54 The pore geometry which might be important is idealized in all cases (including in micro-structures when a morphological image opening depth is used 53 ). Discrepancies between experimental data and PSD model predictions in Figure 4 and 5 might be due to the lack of consideration of the above features in the PSD model.   Catalyst layer.- Figure 6a shows the CL PSD fitting, with the parameters in Table II, against the experimental MIP data. The percentage of hydrophilic pores used in the PSD model is assumed to be 25% since the capillary pressure-saturation relationship for doctor blade CL has not been reported in the literature. A parametric study on hydrophilic volume fraction in CL is provided in a subsequent paper. 46 To validate the ability of the model to reproduce capillary pressuresaturation relationships, the experimental data reported by LaManna et al. 55 was reproduced (see Figure 6c) using the PSD in Figure 6b (provided by LaManna et al. 55 ). Figure 6c also compares the capillary pressure vs. saturation relationships between the doctor blade CL and the CL reported by LaManna et al. 55 The original data reported by LaManna et al. 55 is multiplied by cos(99 • ) in order to account for the assumed contact angle for water in this article. The assumed contact angle of 99 • is based on Ref. 56. The two are remarkably different likely because the CL structures are very different. The pores in the doctor blade CL are centered at 50 nm and 100 nm respectively, with the smaller PSD model assumed to be in hydrophilic phase. Whereas, the majority of pores in the CL reported by LaManna et al. 55 are less than 50 nm. Therefore, the difference in capillary pressure-saturation relationship between the two CLs is due to the different PSDs with higher capillary pressure required to invade the smaller hydrophobic pores in the CL reported by LaManna et al. 55 Figure 6. (a) Experimentally measured PSD using MIP for Doctor blade CL and its fitting, (b) the PSD reported by LaManna et al. 55 57 A detailed model validation with the previous measured and validated PSDs and the CCM fabricated and electrochemically characterized in our laboratory is performed in a subsequent paper. 46 The two operating conditions used in this validation case are 80 • C and a back pressure of 190 kPa, and 40 • C and 150 kPa where the capillary driven flow is dominant. In both cases, the fuel cell is operating at fully humidified condition.
Owejan et al. 57 studied a Gore 5510 MEA with a membrane thickness of 18 μm, and a catalyst loading of 0.4 mg Pt/cm 2 at both anode and cathode. The thickness of the CL is assumed to be 10 μm and the platinum active area is 70 m 2 /g based on values reported previously by Owejan et al. 58 The membrane thermal conductivity for Gore 5510 is reported by Kim et al. 17 to be 0.37 W/(m K). The GDL used is a Mitsubishi Rayon Co. (MRC) 105 with a through-plane thermal conductivity of 0.3 W/m K. The in-plane thermal conductivity is assumed to be ten times higher than through-plane conductivity. The rest of the GDL parameters for example, PSD, porosity and tortuosity, are from the X-ray tomography data reported by Zenyuk et al. 59 Table  III briefly summaries all MPL parameters used in this validation case as well as some of the GDL parameters. The MPL PSD is obtained by subtracting the PSD data of SGL 34BA from SGL 34BC. The other parameters remained the same as mentioned in the Input parameters section. Figure 7 shows that the two-phase model predicted performance is in good agreement with the experimental data at varying operating conditions. At 40 • C where the two-phase flow is dominant, the electrochemical performance with an MPL is significantly better than the one without. The reason for the improved performance is the low MPL thermal conductivity which increases the temperature in the cathode electrode thereby allowing the removal of water in vapor form and minimizing saturation in the CL. The model is capable of capturing the performance change in the two-phase regime with and without MPL. To reduce the number of parameters and sources of uncertainty in the model validation, the authors first focus on studying the fuel cell operation without MPLs in a subsequent paper. 46 The role of MPL will be studied using the two-phase model in future work.

Conclusions
A multi-dimensional, non-isothermal and two-phase membrane electrode assembly model is developed in the open-source fuel cell simulation package OpenFCST. 21 The developed MEA model takes into account gas transport by molecular and Knudsen diffusion, liquid transport by hydraulic permeation and heat pipe effect (phase change), membrane water transport by back-diffusion, thermoosmosis, electro-osmosis drag and heat generation and transport. Double trap and dual path kinetic models are used to predict cathode and anode electrochemical reactions. Local transport losses are accounted for using a local transport resistance as discussed in Refs. 25,26. A dynamic boundary condition at the gas diffusion layer (GDL)/gas channel interface is implemented allowing the model to study both dry and wet conditions as well as the transition between states. The model is used to study fuel cells at varying operating conditions. In order to relate layer micro-structural information to water accumulation, a PSD model, obtained by MIP, is used. The PSD model is used to predict local water saturation based on gas and liquid pressures, and can be used to study the effect of varying pore size and wettability. The micro-structural PSD model is validated by comparing the capillary pressure-saturation and saturation-relative liquid permeability relationships to literature data.
The PSD-based MEA model is validated by comparing numerical predictions against the experimental data in Reference 57. The results show that the two-phase model is able to capture the electrochemical performance variations with and without an MPL at both high and low temperature under fully humidified condition. In the second part of this article, i.e., Ref. 46, the two-phase model is further validated by comparing to: a) experimentally measured electrochemical performance under various operating conditions with CCMs fabricated and characterized in our laboratory; b) liquid water distributions obtained by X-ray in literature; and, c) water fluxes in order to gain a better understanding of the effects of CL and GDL micro-structural changes on fuel cell performance.