A Multi-Fluid Model for Water and Methanol Transport in a Direct Methanol Fuel Cell

: Direct-methanol fuel cell (DMFC) systems are comparatively simple, sometimes just requiring a fuel cartridge and a fuel cell stack with appropriate control devices. The key challenge in these systems is the accurate determination and control of the ﬂow rates and the appropriate mixture of methanol and water, and fundamental understanding can be gained by computational ﬂuid dynamics. In this work, a three-dimensional, steady-state, two-phase, multi-component and non-isothermal DMFC model is presented. The model is based on the Eulerian approach, and it can account for gas and liquid transport in porous media subject to mixed wettability, i.e., the simultaneous presence of hydrophilic and hydrophobic pores. Other phenomena considered are variations in surface tension due to water–methanol mixing and the capillary pressure at the gas diffusion layer–channel interface. Another important aspect of DMFC modeling is the transport of methanol and water across the membrane. In this model, non-equilibrium sorption–desorption, diffusion and electro-osmotic drag of both species are included. The DMFC model is validated against experimental measurements, and it is used to study the interaction between volume porosity of the anode gas diffusion layer and the capillary pressure boundary condition at the anode, and how it affects performance and limiting current density. A.C.O. S.K.K.; funding


Introduction
During the last decades, several types of fuel cells (FCs) have been subject to intensive research, both on component and system-level. While many types of FCs require large volumetric fuel storage or complicated reforming and gas cleaning systems, liquid-fed direct methanol fuel cells (DMFCs) circumvent this by directly oxidizing a fuel consisting of a methanol-water solution. Hence, this type not only offers a high electrical efficiency, it has a high energy density, fast start-up characteristics, and nearly zero recharge time. These characteristics make DMFCs suitable for backup systems and portable power applications in direct competition to batteries. A recent review of the different methanol producing technologies was published by Araya et al. [1].
The system around a direct-methanol fuel cell is very simple, and this is one of the reasons that they have been successfully commercialized already, especially in the enthusiast and leisure market such as camping and sailboats. Figure 1 depicts such an EFOY system by the manufacturer SFC Energy that can provide a maximum charge capacity of 250 Ah/day and a maximum power of 125 W to recharge common types of 12 V or 24 V on-board batteries. According to the manufacturer, a 10 L container filled with pure methanol affords autonomy for up to four weeks, and it can be easily replaced and refilled. While the system of a direct-methanol fuel cell is simple compared to other fuel cell types, the heat and mass transfer processes that occur in these fuel cells are substantially more complex, and this gives rise to the need for detailed computational fluid dynamics analysis. complex, and this gives rise to the need for detailed computational fluid dynamics analysis. A common problem of DMFCs is the high permeation rate of water and methanol through the polymer electrolyte membrane. While methanol crossover results in a loss of fuel and a mixed potential at the cathode [2,3], water crossover causes a flooding of the cathode-backing layer and channels and consequently inhibits oxygen mass transport [4]. The extent of methanol crossover is typically reduced by operating DMFCs with thicker membranes and dilute methanol solutions, whereas the extent of water flooding can be controlled by the air stoichiometry and dew point diagrams [5]. Another two-phase flow phenomenon affecting DMFC operation is gas phase blockage of the anode backing layers and channels. Just as the liquid phase at the cathode can obstruct air mass transport, the gas phase at the anode can obstruct methanol mass transport in the liquid phase [6,7].
An aid in the investigation of these macroscopic transport phenomena governing the electrolyte membrane and backing layers is mathematical modeling. With the ongoing progress in experimental measurements of transport, thermodynamic and electro-chemical properties, the development of detailed mechanistic models of physiochemical processes and increased computational power, incremental improvements are continuously added to mathematical models of DMFCs. In recent years, this has entailed a substantial number of publications in the area of DMFC modeling. In the beginning, these models predominantly covered one-dimensional, multi-component, isothermal and single-phase transport, and were aimed at predicting concentration distributions, methanol crossover and polarization curves [8]. More recently, models have been developed that account for two-phase flow, thermal effects and higher dimensional transport effects as well as detailed membrane transport phenomena. In the following, a short review is given on the progress within these latter areas.
In an early attempt, Wang and Wang [9] investigated the coupling between twophase flow, methanol crossover and the resulting mixed potential in a two-dimensional framework. This was done while ignoring charge transport and thus assuming a uniform overpotential in the anode and cathode catalyst layer (CL) as well as treating the CL as an interface without thickness. Their results underlined the importance of keeping the methanol concentration below 2 M in order to avoid excessive methanol crossover and performance loss. In their work, two-phase flow in the porous media was described using the multiphase mixture formulation, where only one set of governing conservation equations is solved for both thermodynamic phases, as opposed to the two-fluid method where two sets of governing equations are solved. Another feature of the developed DMFC model was the formulation of the capillary pressure in the momentum equations which was described by the dimensionless Leverett J-function [10]. This conveniently expresses the mathematical relationship between capillary pressure, surface tension, viscous permeability, porosity and contact angle of a porous medium in contact with two immiscible fluids as a function of wetting phase saturation. It does so by assuming an idealized hydrophilic or hydrophobic porous medium, consisting of non-connected capillary tubes. One A common problem of DMFCs is the high permeation rate of water and methanol through the polymer electrolyte membrane. While methanol crossover results in a loss of fuel and a mixed potential at the cathode [2,3], water crossover causes a flooding of the cathode-backing layer and channels and consequently inhibits oxygen mass transport [4]. The extent of methanol crossover is typically reduced by operating DMFCs with thicker membranes and dilute methanol solutions, whereas the extent of water flooding can be controlled by the air stoichiometry and dew point diagrams [5]. Another two-phase flow phenomenon affecting DMFC operation is gas phase blockage of the anode backing layers and channels. Just as the liquid phase at the cathode can obstruct air mass transport, the gas phase at the anode can obstruct methanol mass transport in the liquid phase [6,7].
An aid in the investigation of these macroscopic transport phenomena governing the electrolyte membrane and backing layers is mathematical modeling. With the ongoing progress in experimental measurements of transport, thermodynamic and electro-chemical properties, the development of detailed mechanistic models of physiochemical processes and increased computational power, incremental improvements are continuously added to mathematical models of DMFCs. In recent years, this has entailed a substantial number of publications in the area of DMFC modeling. In the beginning, these models predominantly covered one-dimensional, multi-component, isothermal and single-phase transport, and were aimed at predicting concentration distributions, methanol crossover and polarization curves [8]. More recently, models have been developed that account for two-phase flow, thermal effects and higher dimensional transport effects as well as detailed membrane transport phenomena. In the following, a short review is given on the progress within these latter areas.
In an early attempt, Wang and Wang [9] investigated the coupling between twophase flow, methanol crossover and the resulting mixed potential in a two-dimensional framework. This was done while ignoring charge transport and thus assuming a uniform overpotential in the anode and cathode catalyst layer (CL) as well as treating the CL as an interface without thickness. Their results underlined the importance of keeping the methanol concentration below 2 M in order to avoid excessive methanol crossover and performance loss. In their work, two-phase flow in the porous media was described using the multiphase mixture formulation, where only one set of governing conservation equations is solved for both thermodynamic phases, as opposed to the two-fluid method where two sets of governing equations are solved. Another feature of the developed DMFC model was the formulation of the capillary pressure in the momentum equations which was described by the dimensionless Leverett J-function [10]. This conveniently expresses the mathematical relationship between capillary pressure, surface tension, viscous permeability, porosity and contact angle of a porous medium in contact with two immiscible fluids as a function of wetting phase saturation. It does so by assuming an idealized hydrophilic or hydrophobic porous medium, consisting of non-connected capillary tubes. One merit of this approach is the characterization of capillary pressure through easily accessible data. In later work, Berning et al. [11] proposed to correlate the steepness of the Leverett curve with the pore-size distribution, concluding that a wider pore-size distribution of the porous medium results in a steeper capillary pressure curve and is thus desirable. Nonetheless, it should be noted that the validity of the Leverett function is strongly debated due to its limited ability to predict capillary characteristics of the porous media found in FC [12,13].
In the two-dimensional, two-phase and multi-component model by Divisek et al. [14] published around the same time as Wang and Wang [9], a different approach was taken to two-phase modeling. Their model was instead based on the two-fluid approach, which enabled them to account for mixed wettability as well as irreducible saturation in the formulation of the capillary pressure equation. Although a comprehensive model was developed accounting for spatial distributions in the CL, detailed electrochemical reactions, charge transport, dissolved species transport and mixed wettability, the effect of mixed potential was not addressed and validation could not be obtained.
Ge and Liu [15] presented a three-dimensional, two-phase and multi-component liquid-fed DMFC model. The two-phase model was based on the multiphase mixture model, capillary pressure was described by the Leverett J-Function and thermodynamic equilibrium was assumed between phases. Although the model accounted for the spatial resolution of the CL and one-dimensional methanol crossover due to electro-osmotic drag and diffusion, it neglected methanol transport in the CL and thus assumed a uniform distributed parasitic current density across the CL. Their study underlined the improved model predictability of switching from a single-phase to a two-phase flow model at high current densities. In particular, it was shown that the predictability of methanol crossover was improved. Shortly after, Liu and Wang [16] also presented a three-dimensional, isothermal, two-phase model of a liquid feed DMFC, which included proton transport. This model was to some extent a continuation of the work by Wang and Wang [9]. Unlike the previous model, this model included a spatially resolved CL and a non-uniform overpotential across the CL. Compared with Liu and Wang [16], the authors did not account for species transport inside the membrane, but rather formulated source terms at the CL-membrane interfaces assuming a one-dimensional, linear diffusive and convective transport across the membrane. Their model highlighted the effect of the land area on the non-uniform distribution of methanol and current density.
Yang and Zhao [17] presented a two-dimensional, isothermal and two-phase mass transport model for liquid-fed DMFCs. The proposed model accounted for two-phase flow in the porous layers using the two-fluid method and the Leverett J-function. Two different approaches were taken for modeling two-phase flow in the anode and cathode channels. In the cathode, a homogeneous flow model was used, whereas a modified one-dimensional drift-flux model was applied at the anode to account for a difference in phase velocities. The model was later updated by Yang and Zhao [18] to account for species transport and non-equilibrium phase change. This enhancement was used for highlighting the improved predictability by switching from a thermodynamic phase equilibrium condition to a non-equilibrium phase condition while accounting for methanol vapor transport. The developed model was further expanded to a three-dimensional domain by Yang et al. [19]. With this improvement, the effect of the liquid phase saturation at the interface between the channel and the gas diffusion layer was investigated. It was shown that the higher the saturation, the lower was the mass transport resistance for methanol transport. In the work by Xu et al. [20], the issue of non-equilibrium sorption/desorption of water between its dissolved state in the electrolyte phase and a gas-liquid mixture was addressed. This was done using a one-dimensional, isothermal, two-phase model, similar to the mass transport model presented in [17,18,20].
Miao et al. [21] developed a two-dimensional, isothermal, two-phase model of a liquid feed DMFC including electron and proton transport. The employed two-phase model was based on the two-fluid method, the Leverett J-function and non-equilibrium phase change, similar as in [18,19]. Their model, however, accounted for a non-uniform overpotential and parasitic current density, as well as a detailed agglomerate model. The model assumed a fully hydrated membrane and that membrane transport occurred through the liquid phase, similar as in references [9,14]. A paramount issue their model addressed was the treatment of the liquid saturation boundary condition at the GDL-channel interface at the anode. In the majority of DMFC models, a constant liquid saturation interface condition of nearly one is assumed independent of the operation condition [4,16,[18][19][20]22,23]. The underlying basis of this specification has neither been substantiated experimentally nor mechanistically. Miao et al. formulated the GDL-channel interface based on an average liquid saturation between inlet and outlet [21]. Their work highlighted that a significant nonuniform parasitic current density and overpotential distribution can be observed. Later, Miao et al. [22] improved the previously presented model by accounting for homogeneous heat transport, thermal and electrical contact resistance, inhomogeneous compression of the GDL and a more detailed reaction mechanism at the anode.
Yang et al. [24] re-examined the existence of a fixed liquid saturation boundary condition at the GDL/channel interface, which had been used in their previous work [17][18][19][20]23]. It was attempted to remove the GDL/channel interface boundary condition by introducing a liquid water saturation boundary in the CL, and thereby quantifying the liquid distribution in the porous media. The thus obtained boundary condition was currentdensity-dependent.
In a study by Garvin and Meyers [25], a detailed thermodynamic, multicomponent and two-phase model of the anode backing layer was developed. Although only onedimensional, the model combined in detail a non-ideal thermodynamic description of species saturation pressures and a bundle-of-capillaries model of the two-phase flow in porous media. Their model was used for analyzing the promotion of gas removal and a lowering of the methanol concentration in the catalyst layer by tailoring the backing layer properties such as contact angle, permeability, etc., and it was found that the backing layer has to be tailored for specific applications.
He et al. [4] presented a two-dimensional, non-isotherm, two-phase DMFC model including current and proton transport. The developed model is fundamentally similar to the one presented by Miao et al. [21], with two exceptions: the electrolyte model and the anode reaction mechanism. Although the electrolyte model accounts for saturationweighted equilibrium water content in the CL (rather than assuming a fully hydrated membrane as in references [9,14,[17][18][19][21][22][23]), it does neglect a significant difference in the non-equilibrium sorption-desorption kinetic rate of liquids and gases.
In 2013, Bahrami and Faghri [26] published a detailed literature review and highlighted the differences between the various modeling approaches.
Based on the literature review, it can be concluded that within the field of DMFC modeling, the simultaneous presence of hydrophilic and hydrophobic pores has not been taken into account in three-dimensional studies, even though the fraction of hydrophilic pores for SGL-type GDLs ranges between 15-60%, depending on the PTFE content [13,27], and hence should not be neglected. It will be shown below that this is of particular importance for DMFCs where at the anode side the water/methanol mixture reaches the catalyst layer primarily using the hydrophilic pores in the MPL and combination of hydrophilic and hydrophobic pores in the GDL, while the product gas only reaches the channel via the hydrophobic pores in any layer. Hence, in a DMFC model it is of utmost importance to account for this fractional wettability. Moreover, it appears that most modeling attempts for DMFCs that include detailed membrane transport and sorption/desorption in conjunction with detailed two-phase flow in porous media are limited to one-dimensional and sometimes two-dimensional models.
Likewise, most models neglect detailed two-phase flow transport in channels of the anode, even though the amount of gas exiting the anode can be substantial and therefore should not be disregarded. Hence, in the present effort, a boundary condition for the GDL/channel interface based on the channel pressure is prescribed. The anode pressure loss ranges up to several thousand Pascals, and hence constitutes a significant pressure that can force the liquid phase into the hydrophobic pores of the GDL [28]. However, it is doubtful if it is high enough to substantiate a liquid saturation of nearly one at the GDL-channel interface. At the present moment, in order to keep this boundary condition as simple as possible, it is limited to a constant average pressure value for a given channel length. In order to study the importance of a pressure-based boundary condition, two simulation cases with different anode GDL volume porosities are compared and examined in detail with regard to liquid saturation and methanol distribution in the anode. Moreover, focus is put on how these distributions are impacted by fractional wettability of each layer as well as variations in the surface tension.

Mathematical Model
The three-dimensional, steady-state, two-fluid, multicomponent and non-isothermal DMFC model presented here consists of a membrane, anode and cathode electrodes along with adjacent channels. The DMFC model was implemented in the commercial CFD software ANSYS CFX. The two-phase flow model employed is based on the multi-fluid approach, i.e., the model solves one set of transport equations for each phase. The multifluid approach uses the notion of interpenetrating continua, and hence inherently only captures statistical macroscopic phenomena.
Two phases are assumed in the channels, a liquid and a gas phase. In the porous media an additional solid phase is present, and heat transfer in all three phases is included. Moreover, charge transport is added to the porous media. In the membrane, only ion transport, dissolved species and heat transport are modeled.

Computational Grid
The geometry used for this study is depicted in Figure 2. The mesh comprises approximately 52,000 hexahedral elements. The computational domain consists of a membrane electrode assembly (MEA) and two flow channels. Each electrode has a catalyst layer (CL), micro porous layer (MPL) and gas diffusion layer (GDL). The bipolar plates (BP's) are currently neglected. Instead, a temperature and electrical potential boundary condition is applied to channel walls and GDL surface area in contact with BP. ranges up to several thousand Pascals, and hence constitutes a significant pressure that can force the liquid phase into the hydrophobic pores of the GDL [28]. However, it is doubtful if it is high enough to substantiate a liquid saturation of nearly one at the GDL-channel interface. At the present moment, in order to keep this boundary condition as simple as possible, it is limited to a constant average pressure value for a given channel length.
In order to study the importance of a pressure-based boundary condition, two simulation cases with different anode GDL volume porosities are compared and examined in detail with regard to liquid saturation and methanol distribution in the anode. Moreover, focus is put on how these distributions are impacted by fractional wettability of each layer as well as variations in the surface tension.

Mathematical Model
The three-dimensional, steady-state, two-fluid, multicomponent and non-isothermal DMFC model presented here consists of a membrane, anode and cathode electrodes along with adjacent channels. The DMFC model was implemented in the commercial CFD software ANSYS CFX. The two-phase flow model employed is based on the multi-fluid approach, i.e., the model solves one set of transport equations for each phase. The multifluid approach uses the notion of interpenetrating continua, and hence inherently only captures statistical macroscopic phenomena.
Two phases are assumed in the channels, a liquid and a gas phase. In the porous media an additional solid phase is present, and heat transfer in all three phases is included. Moreover, charge transport is added to the porous media. In the membrane, only ion transport, dissolved species and heat transport are modeled.

Computational Grid
The geometry used for this study is depicted in Figure 2. The mesh comprises approximately 52,000 hexahedral elements. The computational domain consists of a membrane electrode assembly (MEA) and two flow channels. Each electrode has a catalyst layer (CL), micro porous layer (MPL) and gas diffusion layer (GDL). The bipolar plates (BP's) are currently neglected. Instead, a temperature and electrical potential boundary condition is applied to channel walls and GDL surface area in contact with BP.

Assumptions
Because of the inherent complexity of the model, a number of simplifying assumptions were made to ensure satisfactory convergence. Care was taken that none of these simplifications had a strong impact on the results discussed below.

•
On the cathode side, the gas phase consists of nitrogen, oxygen and water vapor, whereas the liquid phase only comprises water. The anode gas phase consists of water vapor, methanol and carbon dioxide. The liquid phase is modeled as a binary solution consisting of water and methanol; • Only water and methanol undergo phase change. Species from one phase cannot be dissolved in the other phases and transported; • The cross-over of carbon dioxide through the membrane is currently neglected; • The oxygen reduction reaction (ORR) at the cathode is assumed to produce liquid water. However, when the relative humidity of the adjacent gas phase is below 100%, this water will evaporate immediately; • Gases are assumed to behave ideally, since pressures are low; • In general, all phases share the same pressure field. However, in the CL, MPL and GDL, capillary forces are accounted for that act as so-called body force terms in the momentum equations; • Membrane swelling is currently neglected.

General Governing Equations
The general set of conservation equations for mass, momentum and species are depicted in Equations (1)-(3), respectively.
where the subscripts α and β denote the phases, s the saturation, ε the porosity, ρ the density, . m αβ and . m βα the directionally dependent interfacial mass flows, Y Aα the mass fraction of component A in phase α, τ the stress tensor, p the pressure, and M α the interfacial forces acting on phase α due to the presence of other phases. Moreover, D Aα is the diffusion coefficient of species A in the phase α, and S is the source term. The true phase velocity or intrinsic velocity, U α , is related to the superficial velocity using the expression In addition to mass and momentum, the transport of energy is modeled. However, in the present case, energy transport is simplified by neglecting kinetic energy changes and by solving a homogeneous temperature field. Thus, the energy equation can simply be written as: where H α denotes the static enthalpy, λ α is the thermal conductivity, and N C is the number of species. The molar enthalpy of species i is denoted by h i . The source term S results from all the local overpotentials, both activation and ohmic. These are the general transport equations. Depending on the region of the cell, the general sink and source terms require different expressions which will now be described.

Flow Channels
Inside the channels, momentum is transferred between the phases by an interphase drag force. The current implementation uses the following mixture formulation: where A lg and ρ lg denote interfacial area density and the mixture density, respectively. For two-phase channel flow, the interfacial area density is defined relative to the phase volume fraction and particle size. In this work, a distinction is made between anode and cathode channel due to a difference in the flow morphology. The anode is modeled as a continuous-continuous flow, whereas the cathode is modeled as dispersed-continuous flow. For large gas volume fractions, the gas-liquid flow may approach a continuous-continuous regime, and hence cause a violation of the dispersed gas assumption in the calculation of the interfacial area. Not accounting for this may overpredict interfacial drag in the channel and have a detrimental effect on the convergence behavior. The interfacial area densities of these two types of flows are modeled as follows, where Equation (7) applies to the anode and Equation (8) to the cathode: Here, d is the average particle size, which in our case is assumed to be 20 µm, and s is again the saturation. It should be noted that the abovementioned distinction is only made with regard to interfacial area density.

Porous Media
In the porous media of a fuel cell, the momentum equation is dominated by viscous and capillary forces. One way to describe the effect of viscous forces is through Darcy's law. This equation effectively relates a phase pressure gradient to the phase velocity as follows: where K and k rel denote permeability tensor and relative permeability, respectively, while µ denotes the dynamic viscosity. In general, the relative permeability depends on pore structure, wettability, capillary forces and saturation history [29] and thus is rather complex in nature. However, it is often modeled as being dependent only on saturation. A common approach in FC modeling is a power law correction [9,11,30]: The coefficient q is an empirical constant in the range of 2-5, hence implying a strong interaction with local saturation, and S e is the effective saturation. The effective saturation is a function used to account for the presence of irreducible saturation s irr in a porous medium, i.e., the wetting phase is immobile below a saturation value of s w = s w,irr [31]. The effective saturation S e is defined as follows [32]: With the introduction of Darcy's equation, the presence of capillary forces can be easily accounted for. The capillary pressure is defined as the difference in phase pressure according to: In order to introduce Equations (9), (10) and (13) in the momentum equations, they have to be rewritten as a source terms or so-called body force terms [11,33]. Moreover, if the capillary pressure gradient is known, only one pressure field has to be solved for, which either could be the liquid or gas phase pressure. Assuming that the gas phase pressure is solved for, the resulting momentum source terms are as follows: When modeling the capillary pressure p cap on a macroscopic scale, it is convenient to describe it in terms of the Leverett's function, Equation (16a,b), which combines the fluid properties, material properties and liquid saturation. This approach was first proposed by Udell [34]. However, the Leverett function in its original form can only handle a mixed wettability and no fractional wettability. Hence, in the current study the standard Leverett function is altered by introducing an effective wetting saturation based on the hydrophilic pore fraction, i.e., Equation (17). This approach assumes that the hydrophilic pores are initially filled, and subsequently the hydrophobic pores. Hence, the capillary pressure behaves as follows: p cap < 0 for s < f hi , p cap = 0 for s = f hi , and p cap > 0 for s > f hi" where f hi denotes the hydrophilic pore fraction. A similar implementation has been used in other studies [11,32,35,36]. In references [11,32], it was assumed that only the hydrophobic phase was mobile.
The complete set of equations used to describe the capillary pressure is as follows: where σ denotes surface tension, θ the contact angle, ε the porosity, K the effective permeability and J(S) the Leverett function. In the Leverett function, the term √ K/ε represents the inverse of the characteristic pore radius, 1/r c . This length scale is helpful in evaluating the effect of the capillary pressure.
For methanol-water mixtures, it is important to account for the variation in surface tension due to the methanol molar fraction. The mixture surface tension is described using the following equation [37]: where a and b are empirical coefficients that exhibit a linear dependency on the temperature in the range between 20 • C and 50 • C [37]. In our study, it was assumed that these values are valid at higher temperatures as well, since the temperature dependency is small and linear. The mixture surface tension behaves strongly nonlinear with changes in methanol molar fraction [37], and hence is very important to incorporate. In the present model, both Fickian and Knudsen diffusion are accounted for, similar to the model presented in [38]. The model also accounts for effective properties due to multicomponent diffusion and porous media obstruction. The implemented equations are shown in Table 1.

Parameter Symbol Unit Equation
Gas phase multicomponent correction Temperature and pressure correction of the Fickian diffusivity Porous media correction D e f f Aα where ξ denotes the condensation/evaporation coefficient, and MW is the molecular weight of the species undergoing phase change. Equation (19) exists in various forms. In the following, a simplified version is used where an average temperature is assumed and the liquid phase pressure is given by the saturation pressure. This leads to the following implementation [11,29]: where k xm denotes the convective mass transfer coefficient in [m/s]. The convective mass transfer coefficient was originally derived based on kinetic theory and the assumption of ideal gas behavior; however, the diffusion process needs to be corrected due to the presence of other species by accounting for intermolecular collisions. Further, the ability of the droplet to absorb water molecules needs to be accounted for [32]: where u m denotes the mean molecular speed and Γ the dimensionless uptake coefficient, which accounts for the presence of non-condensable species at the gas-liquid interface.
The interfacial area density depends on the phase change mechanism. In contrast to boiling, evaporation requires a liquid-gas interfacial area in order to occur. This is also in contrast to condensation, which can occur either on a pre-existing liquid layer or on a hydrophilic surface [32]. Inside porous media, it is reasonable to assume that the interfacial area density depends on the pore surface area [32]. Furthermore, the mass transfer rate is corrected for local liquid saturation and porosity: where A pore denotes specific pore surface area, and Γ s the surface accommodation coefficient. Fundamentally, it has to be true in the case of evaporation that A lg → 0 for s → 1 ∨ s → 0 and in the case of condensation that A lg → 0 for s → 1 . Moreover, the surface area exposed to evaporation is only covered by a fraction of a given species. This fraction is assumed equivalent to the species liquid molar fraction. Since methanol and water together make up a non-ideal vapor-liquid mixture, the equilibrium saturation pressure of each species needs to be corrected. At equilibrium, the fugacity of each component in each phase is equal. By introducing an activity coefficient and assuming that the liquid fugacity of component i at standard state is equal to the saturation pressure of its pure component, and assuming an ideal vapor phase (i.e., vapor fugacity coefficient of 1), the vapor pressure of component i can be stated as follows: The activity of any species in the gas or liquid phase can be defined relative to pure standard conditions as: It is important to correct the saturation pressure at low molar fractions of methanol, since it behaves highly non-ideal in this range. The activity coefficients of a vapor-liquid mixture can be estimated via the empirical Non-Random Two Liquid (NRTL) model. Moreover, the saturation pressures of pure water and methanol have a strong non-linear dependence on temperature. To improve numerical robustness, a curve fit to tabulated saturation pressures in an interval between 0 • C and 100 • C is implemented.

Catalyst Layers
The current and ion density distributions are modelled using Ohm's law: where Φ s and Φ m denote the solid potential and membrane potential, respectively, σ the electron conductivity, κ the ion conductivity, and R the charge production or consumption rate. The ion density vector is given as i = −σ∇Φ s and the electric current density vector as j = −κ e f f ∇Φ m . These currents generated in the membrane and solid phase also give reason for ohmic heating, implemented as the source term in the energy equation, Equation (5), where in the regions in question, the source term is S = i 2 /σ and/or S = j 2 /κ. At the interface between the GDL and BP, a voltage drop exists due to contact resistance between the carbon fibers and BP. This contact resistance depends on clamping pressure and PTFE content. The voltage drop is modeled as follows: where R C denotes contact resistance.
The rate of charge production and consumption at the cathode is given based on the oxygen reduction reaction (ORR) and the direct methanol oxidation reaction (DMOR), respectively: In order to determine the volumetric current density at the cathode due the ORR and DMOR, the Butler-Volmer expression is used. This approach assumes a single ratedetermining step.
j e f f 0 is the cathode overpotential and ξ = RT/H 0 is an uptake coefficient based on Henry's law [12], accounting for the difference in the oxygen concentration in the gas and the ionomer phase. Moreover, the exchange current density is corrected for its temperature dependence and active sites being blocked.
The rate of charge production and consumption at the anode is given based on the methanol oxidation reaction: This reaction can be considered to consist of the following elementary reaction steps, as proposed by Gasteiger et al. [39]: The volumetric current density is determined by an expression derived by Meyer and Newman [40] based on the previously shown reaction mechanism, neglecting the reaction in the cathode direction: j e f f 0 where K is an equilibrium constant, η a = Φ s − Φ m − U 0 a is the anode overpotential and C CH 3 OH = (1 − s)C CH 3 OH,l + sC CH 3 OH,g is the corrected methanol concentration. The electrochemical reactions in the anode and cathode CL introduce sink and source terms in the continuity, charge and energy equations given as functions of current density drawn from the DMFC and the number of electrons involved in the consumption or formation of a given species, as listed in Equations (31) and (36).
The simultaneous presence of ORR and DMOR at the cathode electrode gives rise to a mixed potential. By subtracting the virtual parasitic current density from the cathode current density in the charge source terms, the mixed potential is accounted for.
The anode and cathode current densities are defined as follows: where A denotes the cross-section area. The remaining properties used in the electrochemical model are found in Table 2. These values have been assumed and adjusted so that the calculated performance curve provided a good match to experimental results as shown below. It has been assured that these values are within a reasonable range.

Membrane Model
Water and methanol transport in the ionomer is described using dilute solution theory (i.e., Fickian diffusion transport and an electro-osmotic drag term) [19,33,41,42] and including the absorption/desorption term S A : where λ denotes the content of component A in the membrane phase, ρ mem is the dry membrane density (ρ mem = 2000 kg/m 3 ), and EW is the equivalent weight of the membrane (EW = 1.1 kg/mol) [43], while F denotes Faraday's constant (96,485 C/mol). The species content is defined relative to the number of sulfonic acid groups: The species content of Nafion equilibrated with a water-methanol mixture is known to depend non-linearly on species activities and temperature. A few thermodynamic and mechanistic models have been proposed in the literature describing this dependency [40,44,45]. However, none of these models can differentiate between gas and liquid phase equilibration, or what is known as Schroeder's paradox. Essentially, this phenomenon describes how a difference in water content is observed for vapor and liquid at equal activity (i.e., equal chemical potential), even though from a thermodynamic point of view, it should not occur. Some research groups have proposed models explaining Schroeder's paradox based on capillary forces [46] or a difference in surface energy of vapor and liquid water on Nafion [47]. Meanwhile, recent studies have shown the absence of Schroeder's paradox [48,49], if care is put into the thermal history and the equilibration time and activity. However, these observations have not been confirmed by other authors. Thus, in the present modeling effort, two-phase sorption/desorption is accounted for by volume fraction weighting the equilibrium liquid and vapor content.
In the current model, a distinction is made between anode and cathode sorption/desorption. At the cathode, only water exists in both liquid and vapor states; i.e., all methanol that crosses the membrane immediately reacts. Thus, equilibrium sorption of water from liquid and vapor states is given as follows [50]: At the anode, a two-phase mixture of methanol and water exists. For the limiting case of dilute liquid mixtures of water and methanol (i.e., x CH 3 OH < 0.2), simplified relations can be given. Equilibrium sorption of liquid methanol-water solution has been measured by several authors [44,51,52]. Based on these measurements, a constant liquid water sorption is implemented, similar to the one at the cathode. For liquid methanol sorption, a secondorder polynomial fit of the measurements by Ren et al. [52] is used. The methanol content is given as function of methanol activity: λ equi,CH 3 OH,l = 5.635a CH 3 OH,l + 30, 3259a CH 3 OH,l 2 (51) To our best knowledge, no sorption data for gaseous mixtures of methanol and water vapor exist in the literature. Thus, the same dilute behavior of methanol sorption is assumed as in the liquid case, however scaled down relative to the difference in maximum sorption. For water sorption, the equation by Zawodzinski et al. [50] is used, neglecting the effect of methanol on water sorption. λ equi,CH 3 OH,g = 3.585a CH 3 OH,g + 19.2983a CH 3 OH,g 2 (52) where a denotes activity and the subscript equi denotes the equilibrium state. The sorption/desorption source terms are implemented according to Equations (53)- (55). A similar non-equilibrium formulation of the uptake rate has previously been implemented by [12,33,42]. Moreover, in the membrane phase of the cathode CL, only water exists, since methanol becomes oxidized. Hence, sorption/desorption only occurs for water. However, in the anode CL, sorption/desorption can occur for both species: . .
where α denotes area density and k the sorption/desorption kinetic coefficient. The kinetics describing interfacial mass transport of polar vapors and liquids are significantly different in Nafion. Liquid water and methanol sorb approximately ten times faster than their vapor counterpart. Moreover, one needs to differentiate between desorption and sorption kinetics. Desorption is approximately two times faster than sorption for both methanol and water [53][54][55]. The sorption/desorption mechanisms of water and methanol appear to occur similarly [56]. Hence, similar uptake kinetics of methanol and water are assumed. The effective uptake kinetic coefficient is modeled similar to Ge et al. [54]: where V m = EW/ρ m is the partial molar volume of dry membrane, V i = MW i /ρ i,l is the partial volume of species i, and j denotes the phase facing the membrane phase. The diffusivity of water and methanol has during the last decade been subject to intensive discussions due to inconsistencies in the reported values which varied over three orders of magnitude, depending on the measurement technique employed (i.e., mass uptake, permeation and NMR-relaxation) [53]. It was discussed by Majzstrik et al. [53] that these differences were caused by not correctly accounting for membrane swelling and the sorption/desorption phenomenon.
Moreover, local maxima in the Fickian diffusivity as a function of water content have been reported [57,58]. However, it was then shown by our group [59] that the spike in the Fickian diffusivity most likely is a mathematical artifact due to the conversion of the chemical diffusivity into a Fickian diffusivity. We further showed that the diffusivity instead exhibits a transition regime with a sudden increase in diffusivity and afterwards stabilization where a low second-order change with water content is observed, depending on whether membrane swelling is accounted for or not: where λ tp = 2.6225 denotes the transition point and δ ti = 0.8758 the transition interval. This expression reflects the experimental observation by Benziger et al. [60], who used PSGE NMR measurements at long delay times. They showed that at a low water activity or water content, the connection between the hydrophilic pores is low, causing a high tortuosity. Increasing the species content facilitates more interconnected hydrophilic pores with a resulting higher diffusivity. Furthermore, it was shown in [56] that the diffusivity of methanol in Nafion follows the same tendency as water, although diffusing slower than water. The EOD coefficient of water in Nafion has been shown to depict the same discontinuous behavior as the water content. For vapor-equilibrated and pre-dried Nafion, a constant EOD coefficient n d = 1 is observed [61]. Whereas for pre-boiled and liquid-equilibrated Nafion, an EOD coefficient of circa n d = 2.5 has been reported numerous times [61][62][63]. To account for a transition between the two states the following function is used: It was noted in previous work that the application of a constant value for the EOD coefficient leads to a diffusion-only transport mechanism of water inside the membrane [64,65].
In the case that the membrane is equilibrated with a mixture containing water and methanol, it has been shown that both EOD coefficients depend on the molar fraction of methanol present in the bulk mixture [62,63]. At low methanol molar fractions, the EOD coefficient of methanol becomes marginal. It is thought that methanol has a tendency to stick to the polymer backbone structure at low molar fractions. As the methanol molar fraction increases, the total EOD coefficient increases likewise. This phenomenon occurs since the EOD of methanol is more strongly dependent on methanol molar fraction than on the EOD of water. The following equations were fitted to experimental data from [62] for

Boundary Conditions
In order to solve the governing transport equations, a set of boundary conditions is needed for each modeling domain. The used mesh is split up into three sub-domains: anode, membrane and cathode. In addition to interfaces between the sub-domains, extra boundaries are found. All boundaries are marked with a dashed line in Figure 3. The following specifications are given for these boundaries:

Boundary Conditions
In order to solve the governing transport equations, a set of boundary conditions is needed for each modeling domain. The used mesh is split up into three sub-domains: anode, membrane and cathode. In addition to interfaces between the sub-domains, extra boundaries are found. All boundaries are marked with a dashed line in Figure 3. The following specifications are given for these boundaries:  It should be noted that the implemented pressure boundary condition at boundary 1 and 7 implies that the GDL-channel interface capillary pressure is directly specified as a • Boundary 1: This boundary constitutes the interface between the anode channel and anode GDL. It is here that the liquid methanol solution enters and a gas consisting of carbon dioxide, saturated with methanol vapor and water vapor, leaves the GDL. For this boundary, the following is specified: p l,GDL = p l,chan , p g,GDL = p 0 , N e − = 0.
• Boundary 2: This boundary defines the land area in contact with anode GDL, through which electrons and heat is transported: u g = 0, u l = 0, V = V Cell , T = T wall . • Boundary 3: This boundary separates the anode CL from the membrane and is impermeable to gas and electrons: N e − = 0, N g,CH 3 OH = 0, N g,H 2 O = 0, N g,CO 2 = 0. • Boundaries 4 and 5: Both boundaries are treated as mirror planes, meaning that any change of any variable in the x-direction has to be equal to zero: dϕ/dx = 0. • Boundary 6: Similar to boundary 3, this boundary separates the CL from the membrane at the cathode, and is impermeable to gas and electrons: N e − = 0, N g,N 2 = 0, N g,H 2 O = 0, N g,CO 2 = 0, N g,O 2 = 0. • Boundary 7: At this boundary, the liquid phase exits the GDL and enters the channel as droplets. Simultaneously, air enters the GDL from the channel. For this boundary, the following is specified: p l,GDL = p l,chan , p g,GDL = p 0 , N e − = 0. • Boundary 8: This boundary defines the cathode land area in contact with GDL, through which electrons and heat is transported: It should be noted that the implemented pressure boundary condition at boundary 1 and 7 implies that the GDL-channel interface capillary pressure is directly specified as a pressure difference between the liquid phase channel and GDL gas phase pressure, as follows: This approach is in contrast to earlier attempts where either a constant liquid saturation or a channel-averaged liquid saturation has been specified. The reason for altering this boundary condition is to obtain a more physically correct boundary condition that is based on pressure rather than saturation. Interestingly, what this boundary condition effectively enforces is that the liquid phase becomes pushed into the GDL due to an overpressure in the channel relative to the gas phase pressure inside the GDL.

Results and Discussion
In this section, the mathematical model is validated and detailed results are presented, including the distributions of liquid methanol concentration, liquid saturation, fluid temperature and dissolved methanol and water content in the electrolyte phase.
In the one-dimensional variable distributions, attention is on the impact of the anode capillary pressure boundary condition at the channel-GDL interface. In order to investigate the importance of this boundary condition, two cases are depicted in each distribution. Each case represents a specific porosity of the anode GDL. The porosity affects both the relative permeability as well as the characteristic pore radius used in the capillary pressure function [66].
In Case 1, a porosity of 0.75 is specified and in Case 2, a porosity of 0.8. This is equivalent to a permeability of 1.2 × 10 −11 m 2 and 3.2 × 10 −11 m 2 , and a characteristic pore size of 8 µm and 10 µm, respectively.
In the three-dimensional plots, emphasis is solely on discussing some of the phenomena the model is able to capture and how they affect some variable distributions. The various input parameters for this model are listed in Tables 3 and 4. Moreover, the electrical contact resistance between the GDL and the bipolar plate was 20 mΩ-cm 2 [67], and the thermal contact resistances were 1.5 × 10 −4 K-m 2 /W [68].

Predicted Polarization Curve
The presented model with the base case parameters has been validated against experimental measurements of the voltage to current density relation. Figure 4 shows the comparison for a bipolar plate temperature of 348 K (i.e., 75 • C at the bipolar plates, see Table 3) The electrochemical parameters in Table 2 were curve-fitted to match these experimental results.

Predicted Polarization Curve
The presented model with the base case parameters has been validated against experimental measurements of the voltage to current density relation. Figure 4 shows the comparison for a bipolar plate temperature of 348 K (i.e., 75 °C at the bipolar plates, see Table 3) The electrochemical parameters in Table 2 were curve-fitted to match these experimental results. With this, the predicted performance is in excellent agreement with the experimental data. At low current densities, in the activation overpotential region, the model captures the change in voltage, especially in the important low current density regime that is dominated by cross-over of methanol. This initial change in voltage is dominated by a change in anode overpotential. The cathode overpotential loss is already high due to methanol crossover. As the current density increases, the change in voltage becomes linear, indicating the beginning of the ohmic region. It seems that the predicted curve falls steeper than the experimentally measured one. This suggests a slightly higher resistance due to either contact resistance, ionic or electron transport. However, it could also be due to an underestimation of the anode overpotential loss. Further improvements may be obtained by slightly adjusting these abovementioned parameters, but overall, the agreement is deemed very satisfactory. With this, the predicted performance is in excellent agreement with the experimental data. At low current densities, in the activation overpotential region, the model captures the change in voltage, especially in the important low current density regime that is dominated by cross-over of methanol. This initial change in voltage is dominated by a change in anode overpotential. The cathode overpotential loss is already high due to methanol crossover. As the current density increases, the change in voltage becomes linear, indicating the beginning of the ohmic region. It seems that the predicted curve falls steeper than the experimentally measured one. This suggests a slightly higher resistance due to either contact resistance, ionic or electron transport. However, it could also be due to an underestimation of the anode overpotential loss. Further improvements may be obtained by slightly adjusting these abovementioned parameters, but overall, the agreement is deemed very satisfactory.

Comparison of Two Cases
In the following, the two cases with different GDL porosities and different characteristic pore radii as mentioned above are examined in detail. The distribution of liquid saturation in the anode electrode is shown in Figure 5a. It is evident that the two cases nearly coincide in the CL and MPL, whereas a clear distinction is possible in the GDL. Even though a large difference in liquid saturation level appears, the same trend of a gas build-up under the land area is visible. This build-up of gas is not only a reflection of transport resistance, but also due to a difference in surface tension. Under the land, the molar fraction of methanol is lower than under the channel, as seen in Figure 5b. This difference in molar fraction imposes a difference in surface tension under land as this is highly dependent on methanol molar fraction. Moreover, this results in a difference in liquid saturation in order to balance the capillary pressure. Another phenomenon that the model can capture is the saturation jump condition at the CL-MPL and MPL-GDL interface. This jump condition in liquid saturation is caused by an abrupt change in porous media properties such as porosity, permeability and hydrophilic pore fraction. By having a highly hydrophobic MPL with small pores, the liquid saturation level can be kept fairly low. This effectively decreases methanol diffusivity thereby minimizing the methanol crossover rate. Unfortunately, this also affects the extent of the anode overpotential loss when the limiting current sets in.
The use of micro-pores has another interesting benefit. In the CL and MPL, this pushes the liquid saturation close to the hydrophilic pore fraction. By having a small characteristic pore size of 0.163 μm and 0.426 μm, respectively, it takes a large pressure for the liquid phase to overcome the repelling capillary forces and intrude the hydrophobic pores. Meanwhile, the characteristic pore size of the GDL is much larger. Consequently the pressure requirement for intruding these hydrophobic pores is less. As it appears from Figure 5, by decreasing the porosity of the GDL slightly, the saturation level is pushed considerably closer to the hydrophilic pore fraction. Again, this occurs due to the increasing capillary forces, which the liquid phase has to overcome. It is evident from these two cases that a too small characteristic pore size of the GDL significantly lowers the liquid saturation level. In turn, this leads to a decrease in the diffusivity of methanol in the liquid Another phenomenon that the model can capture is the saturation jump condition at the CL-MPL and MPL-GDL interface. This jump condition in liquid saturation is caused by an abrupt change in porous media properties such as porosity, permeability and hydrophilic pore fraction. By having a highly hydrophobic MPL with small pores, the liquid saturation level can be kept fairly low. This effectively decreases methanol diffusivity, thereby minimizing the methanol crossover rate. Unfortunately, this also affects the extent of the anode overpotential loss when the limiting current sets in.
The use of micro-pores has another interesting benefit. In the CL and MPL, this pushes the liquid saturation close to the hydrophilic pore fraction. By having a small characteristic pore size of 0.163 µm and 0.426 µm, respectively, it takes a large pressure for the liquid phase to overcome the repelling capillary forces and intrude the hydrophobic pores. Meanwhile, the characteristic pore size of the GDL is much larger. Consequently, the pressure requirement for intruding these hydrophobic pores is less. As it appears from Figure 5, by decreasing the porosity of the GDL slightly, the saturation level is pushed considerably closer to the hydrophilic pore fraction. Again, this occurs due to the increasing capillary forces, which the liquid phase has to overcome. It is evident from these two cases that a too small characteristic pore size of the GDL significantly lowers the liquid saturation level. In turn, this leads to a decrease in the diffusivity of methanol in the liquid phase, hence limiting current density and performance. For Case I and Case II, this gives a current density of 0.17 A/cm 2 and 0.28 A/cm 2 , respectively, a surprisingly large difference for such a comparatively small change in the physical properties.
A clear distinction between the two cases can also be seen in the liquid phase methanol distributions, shown in Figure 5b. As discussed above, a lower liquid saturation level in the GDL causes a decreased methanol transport rate. Hence, a much steeper methanol concentration gradient is visible in the through-plane direction of the GDL for Case I. The same applies to the in-plane direction towards the middle of the land. However, not only a decrease in saturation plays a role here. A decrease in porosity and a corresponding increase in tortuosity decrease the effective diffusivity of methanol. Meanwhile, it should be noted that the distribution of methanol in the liquid phase is furthermore dependent on the rate of methanol evaporation and condensation. Even though the net transport of gas phase is in the direction of the channel, that of methanol vapor may move toward the CL and help improve performance or the limiting current density. However, the extent of this effect on performance is fairly low, as the mass transport rate in the gas phase is lower compared to that of the liquid phase.
The fluid temperature distribution in the anode electrode is shown in Figure 6a. Again, a clear difference can be seen between the two cases. In Case II, where the higher current density is produced, the temperature is also highest; this difference is due to the rate of heat production associated with catalytic burning of crossover methanol, the ORR and MOR reactions, and electron and ion transport. Importantly the methanol crossover generates a high production of heat and temperature rise. This difference in the temperature distribution is also interesting since it affects the sorption/desorption in the electrolyte phase of the CL, species crossover, and phase change. The fluid temperature distribution in the anode electrode is shown in Figure 6a. Again, a clear difference can be seen between the two cases. In Case II, where the higher current density is produced, the temperature is also highest; this difference is due to the rate of heat production associated with catalytic burning of crossover methanol, the ORR and MOR reactions, and electron and ion transport. Importantly the methanol crossover generates a high production of heat and temperature rise. This difference in the temperature distribution is also interesting since it affects the sorption/desorption in the electrolyte phase of the CL, species crossover, and phase change. The decrease in performance between cases with a different porosity is visible from cross sectional area-averaged gas volume fraction in the channel, depicted in Figure 6b. For Case II, a much larger amount of gas is predicted in the channel. Inherently, this occurs since a higher current density means a higher production rate of carbon dioxide. However, the exact amount of gas leaving the GDL is a complex matter to predict, since it depends on the rate of methanol and water evaporation, which in turn depends on temperature.
The three-dimensional distribution of the gas volume fraction in the anode gas diffusion electrode and channel is shown in Figure 7a. A fairly even distribution is seen inside the electrode along the channel length. This occurs in opposition to the channel The decrease in performance between cases with a different porosity is visible from cross sectional area-averaged gas volume fraction in the channel, depicted in Figure 6b. For Case II, a much larger amount of gas is predicted in the channel. Inherently, this occurs since a higher current density means a higher production rate of carbon dioxide. However, the exact amount of gas leaving the GDL is a complex matter to predict, since it depends on the rate of methanol and water evaporation, which in turn depends on temperature.
The three-dimensional distribution of the gas volume fraction in the anode gas diffusion electrode and channel is shown in Figure 7a. A fairly even distribution is seen inside the electrode along the channel length. This occurs in opposition to the channel height, where a gradual increase is visible. This increase was also observed in Figure 5, and it occurs because gas continuously enters the channel from the electrode. When taking a closer look at the middle of the channel, it can be observed that the gas phase is pushed away from the GDL-channel interface and is concentrated at a short distance from the GDL-channel interface. This trend is quite different near the channel corner facing the GDL. Here, a large gas pocket is forming, which gradually increases along the channel length, obstructing more and more liquid phase transport. This phenomenon occurs due to the velocity distribution around the corner rather than surface tension forces, as they are not included at the present stage. In Figure 8, the methanol and water content distributions of the electrolyte phase are shown. In both figures, a steeper gradient can be seen in the CL, caused by a decrease in ionomer content in comparison to the membrane. Decreasing the ionomer content in the CL effectively lowers the diffusivity of methanol and, unfortunately, also the ion conductivity. However, this gradient is not only caused by a decrease in the ionomer content, it is partially also due to EOD and non-equilibrium uptake that balances diffusion. The EOD is quite large in the CL, since the current density changes in the through-plane direction, imposing a large potential field gradient. In Figure 7b, an overview is given of the methanol concentration distribution along the channel for Case 1. Again, a fairly even distribution of methanol is seen along the channel, which results from the high stoichiometry. It is only close to the GDL-channel interface that a gradient in the methanol concentration is visible. This thinning out of the methanol concentration inside the catalyst layer causes a high concentration overpotential, but it also reduces the methanol crossover effect, which is one of the biggest loss mechanisms in DMFCs.
In Figure 8, the methanol and water content distributions of the electrolyte phase are shown. In both figures, a steeper gradient can be seen in the CL, caused by a decrease in ionomer content in comparison to the membrane. Decreasing the ionomer content in the CL effectively lowers the diffusivity of methanol and, unfortunately, also the ion conductivity. However, this gradient is not only caused by a decrease in the ionomer content, it is partially also due to EOD and non-equilibrium uptake that balances diffusion. The EOD is quite large in the CL, since the current density changes in the through-plane direction, imposing a large potential field gradient.
shown. In both figures, a steeper gradient can be seen in the CL, caused by a decrease in ionomer content in comparison to the membrane. Decreasing the ionomer content in the CL effectively lowers the diffusivity of methanol and, unfortunately, also the ion conductivity. However, this gradient is not only caused by a decrease in the ionomer content, it is partially also due to EOD and non-equilibrium uptake that balances diffusion. The EOD is quite large in the CL, since the current density changes in the through-plane direction, imposing a large potential field gradient. Inherently, the distributions of methanol and water are rather different in the membrane. As seen from Figure 7a, the methanol content approaches zero near the interface between the cathode CL and membrane, whereas the distribution of water content is again fairly even. The methanol content has to approach zero at the interface, since methanol is catalytically burned when getting in contact with air. The water content level is dependent Inherently, the distributions of methanol and water are rather different in the membrane. As seen from Figure 7a, the methanol content approaches zero near the interface between the cathode CL and membrane, whereas the distribution of water content is again fairly even. The methanol content has to approach zero at the interface, since methanol is catalytically burned when getting in contact with air. The water content level is dependent on the air's relative humidity of the cathode, which is a function of water crossover, the water production rate due to ORR and DMOR, as well as air stoichiometry.
Finally, it should be noted that the methanol and water content are highly dependent on the hydrophilic pore fraction of the anode CL. Since the liquid saturation is nearly identical to the hydrophilic pore fraction, and the equilibrium methanol and water content are highly dependent on the liquid saturation, and since the liquid phase has a higher equilibrium content and sorption rate, the hydrophilic pore fraction effectively determines the anode water and methanol content.

Conclusions
Direct methanol fuel cell systems are typically simple in comparison with other types of fuel cells. However, the underlying heat and mass transfer mechanisms are exceedingly complex, and computational fluid dynamics modelling can help to shed light into how to improve system performance. In this work, a steady-state, three-dimensional, twofluid, multi-component and non-isothermal liquid-fed direct methanol fuel cell model was presented and compared against experimental measurements in good agreement. The developed model was used for investigating the interaction between volume porosity of the anode GDL and the capillary pressure boundary condition on fuel cell performance. The simulation results indicate that while keeping the force of liquid phase intrusion into the GDL constant, the change in GDL porosity plays a significant role in performance. This would not have been observed for a fixed saturation condition as it implicitly adjusts the force required for intruding hydrophobic pores when changing the characteristic pore size radius of the GDL and hence its capillary pressure. It was further shown that accounting for variations in surface tension due to methanol decreases the prediction of liquid saturation under the land as opposed to under the channel, which pronounces the resistance to methanol diffusion under the land area. Moreover, the importance of accounting for the hydrophilic pore fraction was underlined. The smaller the characteristic pore size becomes for highly flooded porous media, the closer the liquid saturation is pushed to the hydrophilic pore fraction, making methanol transport losses more prominent.