## Abstract

An anomaly in the modelling of two-phase flow in the porous cathode gas diffusion layer (GDL) of a polymer electrolyte fuel cell is investigated asymptotically and numerically. Although not commented on previously in literature, the generalized Darcy model used most commonly leads to the surprising prediction that a hydrophilic GDL can lead to better cell performance, in terms of current density, than a hydrophobic one. By analysing a reduced one-dimensional steady-state model and identifying the capillary number as a small dimensionless parameter, we find a potential flaw in the original model, associated with the constitutive relation linking the capillary pressure and the pressures of the wetting and non-wetting phases. Correcting this, we find that, whereas a hydrophilic GDL can sustain a two-phase (gas/liquid) region near the water-producing catalytic layer and gas phase only region further away, a hydrophobic GDL cannot; furthermore, hydrophobic GDLs are found to lead to better cell performance than hydrophilic GDLs, as is indeed experimentally the case.

## 1. Introduction

The past decade has seen a proliferation of modelling activity on the polymer electrolyte fuel cell (PEFC); an important subset of this activity is the modelling of the two-phase flow that occurs in the gas diffusion layer (GDL) on the cathode [1,2]. The prevailing approach uses a generalized form of Darcy's law, which has been widely used over the past several decades to analyse the movement of oil and water in soils and porous rock [3]. Applied to water transport in PEFCs, the Darcy model characterizes the response of the porous material by the capillary pressure, the gas and liquid phase relative permeabilities, and the effective gas diffusion coefficient, all of which depend on the fraction of the local pore volume occupied by liquid water. Although there are other methods for modelling two-phase flow in the GDL—for example, lattice Boltzmann simulation [4]—the generalized Darcy model is still frequently used [5] and therefore merits closer scrutiny.

In the context of PEFCs, it was first used by Wang *et al.* [6], and then subsequently by Nam & Kaviany [7] and Pasaogullari & Wang [8]. At the heart of the model is the capillary pressure, *p*_{c}, which is a function of the water saturation, *s*; a commonly used expression is
1.1where *σ* is the surface tension, *θ*_{c} is the wetting angle, *γ* is the porosity, *κ* is the absolute permeability and *J* is the Leverett function. From (1.1), which previously had been used for hydrophilic porous rocks and soils for which 0^{°}≤*θ*_{c}<90^{°}, Pasaogullari & Wang [8] postulated the use of
1.2as an extension for hydrophobic media. Although there are alternative forms for *p*_{c} [9,10], equation (1.1) appears to be the only one that contains the wetting angle explicitly; consequently, allowing it to vary in a model may give an indication of how the degree of hydrophobicity, or polytetrafluoroethylene (PTFE) content, of a GDL affects, for example, the performance of the cell. Experimentally, it is already well established that increased PTFE content can have a beneficial effect on cell performance [11–15], although too much PTFE loading can lead to flooding at the catalyst layer [16]. As regards modelling, although there is some evidence that a GDL treated with PTFE does not have uniform wettability [17], but will have both hydrophilic and hydrophobic pores, it remains nevertheless common practice to assume that *θ*_{c} is constant throughout the GDL. A related issue is whether *θ*_{c} should even be present at all in equation (1.2): Slattery [18] indicated that it should, whereas Litster & Djilali [19] pointed out that Anderson [20] had previously indicated that the capillary pressure ought not to be related to the wetting angle in the manner given in (1.2), as have Gostick *et al.* [21]. However, it is clear that omitting *θ*_{c} completely would make it difficult to judge the effect of hydrophobicity in a model; the only option would be to use experimental data for *p*_{c} for two different GDLs made of the same material, typically carbon cloth or paper, but with differing PTFE content.

Furthermore, it has never been conclusively established what effect *θ*_{c} has on model results for the GDL; for the most part, this has been because, although often quoted, equation (1.2) is only ever used for one value of *θ*_{c} in any one investigation. Indeed, we have found only three exceptions [8,22,23]. Senn & Poulikakos [22] compared the water saturation at the catalyst layer for different values of *θ*_{c}>90^{°}, and concluded that a more hydrophobic GDL will lead to lower water saturation at the catalytic layer, and hence better cell performance, because the reaction sites are less blocked; similar conclusions are drawn by Pasaogullari & Wang [8] and Song *et al.* [23], although these articles also contain results for *θ*_{c}<90^{°}. However, from a closer comparative inspection of figs 8 and 9 in Pasaogullari & Wang [8] and figs 14a,c in Song *et al.* [23], it is evident that the saturation at the catalyst/GDL interface decreases as *θ*_{c} decreases: in the former, the saturation at the GDL/catalyst interface is greater when *θ*_{c}=91^{°} than when *θ*_{c}=81^{°}, whereas in the latter, its steady-state value is greater when *θ*_{c}=100^{°} than when *θ*_{c}=0^{°}. Ultimately, this would mean that the current density should increase for more hydrophilic GDLs, which appear to go against conventional wisdom that hydrophilic GDLs are not as able to transport liquid water away from the catalytic reaction layer, to the detriment of cell performance. In general, the absence of model results for *θ*_{c}<90^{°} is rather conspicuous, especially because the theory for *θ*_{c}>90^{°} for PEFC GDLs grew out of previously well-established theory for hydrophilic porous media; thus, the purpose of this paper is to address this issue.

In this paper, we analyse in more detail the simplest possible case, by considering a one-dimensional steady-state isothermal model for the GDL. Even a one-dimensional model is not without practical meaning, because it corresponds to a flow with high stoichiometry; moreover, it forms the basis for two- and three-dimensional models, transient models and models that include a hydrophobic microporous layer that is often added adjacent to the active layer. The paper is organized as follows. The model equations are presented §2, and then non-dimensionalized in §3. A reduced model is analysed asymptotically in §4, and these results are complemented with computations in §5, wherein comparison is also made with earlier work. Conclusions are drawn in §6.

## 2. Model equations

We use the version of multi-fluid model formulation given by Vynnycky [24]. There is essentially no difference between this version of the multi-fluid model and that used by others [25–27], as was shown numerically by Vynnycky [24]; it does, however, remove the more *ad hoc* nature of the way that the inter-phase transfer is treated, thereby making the analysis mathematically more transparent. We consider an isothermal GDL of thickness *H*, in some part or all of which there is two-phase flow, as shown in the electronic supplementary material, figure S1.^{1} The electrochemical reaction
2.1occurs at the catalytic layer at *y*=0, whereas *y*=*H* denotes the GDL/channel interface. Although perhaps not obvious from the outset, there turn out to be three cases to consider:

I. two-phase flow for 0≤

*y*≤*λ*and one-phase flow for*λ*≤*y*≤*H*;II. two-phase flow for 0≤

*y*≤*H*; andIII. one-phase flow for 0≤

*y*≤*λ*and two-phase flow for*λ*≤*y*≤*H*.

In cases I and III, *λ* is not known *a* *priori*, and must be determined as part of the solution; in case II, *λ*=1 implicitly. Case I corresponds to the situation where the liquid water that is produced in reaction (2.1) does not occupy the entire GDL. Case II arises when liquid water is produced in such large quantities that it occupies the entire GDL; this may occur in tandem with the presence of liquid water at the GDL/channel interface. Case III arises when liquid water is present at the GDL/channel interface; this case has received less attention earlier than the others, but we will explain later how it arises. In addition, and as in all earlier papers on two-phase flow in the GDL of a PEFC, we assume negligible quantities of oxygen and nitrogen in the liquid phase.

### (a) Case I

Starting with steady-state conservation equations for both gas and liquid phases for 0≤*y*≤*λ*, we have
2.2and
2.3where denotes the mass flux for species *i*, *ρ*^{(l)} is the liquid density, *v*^{(l)} is the liquid velocity and is the interface mass transfer of water between the gas and liquid phase. Adding the three equations in (2.2) to (2.3) gives
2.4where *ρ*^{(g)} is the gas mixture density, and *v*^{(g)} is the mass-averaged velocity of the gas phase; this eliminates so that only the first two equations in (2.2) and (2.4) need be considered, although can be computed *a posteriori* if necessary [24]. In turn, we have
2.5and
2.6in (2.5), *κ*^{(g)}_{rel}(*s*) is the gas relative permeability, *κ*^{(l)}_{rel}(*s*) is the liquid relative permeability, *μ*^{(g)} and *μ*^{(l)} are, respectively, the gas and liquid phase dynamic viscosities, is the mass fraction and describes the diffusion-driven transport. For the latter, we use
2.7where is the (*i*,*j*)-component of the multi-component Fick diffusivity matrix—which is modified by (1−*s*)^{ϖ} to take account of the liquid phase, where typically *ϖ*=1—and denoting the mole fraction of species *j*, is related to by
in turn, this introduces the relative molecular weights of water, nitrogen and oxygen (*M*_{H2O},*M*_{N2} and *M*_{O2}, respectively) and the mixture molecular weight, *M*^{(g)}, given by
2.8For a ternary system, these are related to the multi-component Maxwell–Stefan diffusivities (*D*_{ij}) through, for *i*,*j*,*k*=H_{2}O, N_{2}, O_{2},
2.9and
2.10Thence, for gases at low density, the multi-component Maxwell–Stefan diffusivities, *D*_{ij}, can be replaced with the binary diffusivities, , for all pairs of species in the mixture; explicit expressions for these, based on the Chapman–Enskog theory, can be found in Bird *et al.* [28]. In addition to the above, we must have
2.11there is no differential equation for because it takes the saturation value, i.e.
2.12where
2.13with *T* as the temperature and *t*=*T*−273.15.

For *λ*≤*y*≤*H*, where *s*≡0, we have
2.14
2.15
and
2.16as well as equation (2.11).

Constitutive relations are required for *ρ*^{(g)}, *κ*^{(g)}_{rel} and *κ*^{(l)}_{rel}. For *ρ*^{(g)}, we use the ideal gas law,
2.17where is the universal gas constant. We do not yet take any particular forms for *κ*^{(g)}_{rel}(*s*) and *κ*^{(l)}_{rel}(*s*), other than to require that
A further relation is required to relate *p*^{(g)} and *p*^{(l)}; this is done via the capillary pressure, *p*_{c}, which is itself a function of the saturation. By definition,
2.18where *p*_{nw} and *p*_{w} denote the pressures of the non-wetting and wetting phases, respectively [3]; so,
2.19In fact, our use of equation (2.19), rather than
2.20constitutes the decisive difference, when compared with earlier models.

The boundary conditions are at *y*=0,
2.21and
2.22where *α* is the number of water molecules dragged through the membrane by each proton, *I* is the current density and *F* is Faraday's constant (96487 C mol^{−1}); at *y*=*H*,
2.23In particular, *I* is normally chosen to be of the form
2.24where *I*_{0} is current density that one would expect in the absence of liquid water. For this, we take
2.25where *a* is the exchange current density, and *η* is the cathodic overpotential; equation (2.25) was used by Vynnycky [24], and was based on that in Natarajan & van Nguyen [29].

Following Vynnycky [24], we have at *y*=*λ*, i.e. the interface between one- and two-phase flow:

— zero liquid saturation and the gaseous water mole fraction in the one-phase region equal to the saturation value: 2.26

— continuity of gas pressure, gaseous oxygen mole fraction, oxygen gas flux and nitrogen gas flux: 2.27

— continuity of water flux summed over all phases, which reduces to and then 2.28

in equations (2.27) and (2.28), denotes the difference in the value of a function above and below *y*=*λ*.

## 3. Non-dimensionalization

We write
where [*p*_{c}]=*σ*/*κ*^{1/2}, and [*v*^{(l)}],[*v*^{(g)}] and Δ*P* are to be determined. It turns out to be appropriate to take
which, in turn, requires a scale for [*v*^{(g)}]. There are two possibilities for [*v*^{(g)}],
The first of these is appropriate at lower current densities, whereas the second is appropriate at high or limiting current density when the electrochemical reaction is diffusion-limited; it turns out to be simpler to use the second of these. For [*I*], we take
3.1whereas [*η*] varies between 0 and 0.6 V.

### (a) Governing equations

#### (i) Case I

For 0≤*Y* ≤*Λ*, where *Λ*=*λ*/*H*,
3.2and
3.3where
3.4
3.5and *Π*=Δ*P*/*p*^{out}. Also,
3.6and
3.7where with
3.8where *Ca* is the capillary number. Combining (3.2), the second equation in (3.6) and (3.7) eliminates *V* ^{(l)} and *P*^{(l)} to give
3.9instead of (3.2); in (3.9), the prime denotes differentiation with respect to *s*, whereas the ± sign refers to the case respectively.

For the region of one-phase flow, *Λ*≤*Y* ≤1, we have just
3.10
3.11
and
3.12

### (b) Boundary conditions

#### (iv) Case I

At *Y* =1,
3.13

At *Y* =0,
3.14
3.15
and
3.16where *Ω*=[*I*][*M*^{(g)}]/*F*[*ρ*^{(g)}][*v*^{(g)}].

At *Y* =*Λ*,
3.17
3.18
and
3.19where .

### (c) Dimensionless parameters

There remain six dimensionless parameters: *Π*,*χ*,*s*^{in}. With
we find that
leading to
with V, we have, in addition, that . These values suggest a reduced model for the regime . Although none of the other remaining parameters are particularly small, we will proceed under the assumption that *Π*≪1, as this enables us to demonstrate the qualitative behaviour of the model that we later show numerically to be valid even when this assumption is not made.

## 4. Reduced model

Proceeding with , *Π*≪1, there are several sub-cases to consider, depending on whether or not *s*^{in}=0, and on whether which corresponds to *θ*_{c}≈ 90^{°}. We consider each of these cases in turn.

Before writing down the reduced model equations, it is worth pointing out the fact that could lead us to suppose that (3.9) can be simplified to give, at leading order,
4.1indeed, this is feasible in the case when *s*^{in}>0, as is shown in §4*a*(ii). However, if this is not the case, meaning that either there is only two-phase flow, so that *s*=0 at *Y* =1, or there is both one- and two-phase flow, so that *s*=0 at *Y* =*Λ*(<1), then the fact that *κ*^{(l)}_{rel}(0)=0 indicates that this would not be correct. Instead, it is better to define a new variable,
4.2which we then rescale according to where we suppose that *Σ*∼*O*(1); furthermore, if *P*^{′}_{c}(0)=0, it is necessary to use
4.3although we will assume henceforth that Equation (4.2) was used by Vynnycky [24], where it was found to be of use in treating the numerical singularity that arises precisely because *κ*^{(l)}_{rel}(0)=0. For example, we see that if *κ*^{(l)}_{rel}(*s*)∼*s*^{n} as where *n*>0, we will have It is noteworthy that the value of *n* determines how large *s* will become; for example, for *n*=3, which is the largest value that is usually ever taken, To proceed further, we suppose that which further suggests that we can write
to leading order in This leads to a decoupling of the equation system, as shown below.

With *Π*≪1, we have [*ρ*^{(g)}]=[*M*^{(g)}]*p*^{out}/*RT*, giving
4.4in addition, is constant. For 0≤*Y* ≤*Λ*, the governing equations are
4.5and
4.6whereas for *Λ*≤*Y* ≤1,
4.7and
4.8As for the boundary conditions, the equations in (3.13) remain unchanged, but (3.14)–(3.16) become
4.9
4.10
and
4.11Hence, *Σ* can be found *a posteriori*, i.e. after *V* ^{(g)} and have been found. From (4.10), we obtain, for 0≤*Y* ≤*Λ*,
4.12substituting into (4.9) gives
4.13where
4.14with thence,
4.15where *C* is a constant of integration. Then, integrating (4.11) once and using *Σ*(*Λ*)=0 gives
4.16

In summary, we have found the analytical form of the leading-order solution in the two-phase region up to three constants—*C*,*I*_{0}, *Λ*—and we would then need to find the solution in the one-phase region. However, we can infer the behaviour of the solution as a function of *θ*_{c}. First, from equation (4.16), we have
4.17of apparent importance is the sign of the right-hand side in (4.17). According to the definition for *p*_{c} in equation (1.2), *P*^{′}_{c}(0)<0 for all *θ*_{c}. Furthermore, using (2.9) and (2.10), and because and the integral in (4.17) will be strictly positive. In addition, because literature values of *α* are typically positive, varying from 0.3 to 1.7 [30], the first term in the numerator of (4.17) will also be positive. Hence, *Σ*(0)>0 for *θ*_{c}<90^{°} and *Σ*(0)<0 for *θ*_{c}>90^{°}. Moreover, because *Σ*(0)<0 would be physically unrealistic, the interpretation is that this would not be an admissible solution.

Next, we consider the case when *s*^{in}>0.

First, we set to obtain, at leading order in ,
4.18and
4.19subject to, at *Y* =0,
4.20
4.21
and
4.22This time,
4.23leading to
4.24and thence
4.25Thus, the behaviour of *S* here will be similar to that of *Σ* in §4(*a*), in the sense that *S*(0)>0 for *θ*_{c}<90^{°} and *S*(0)<0 for *θ*_{c}>90^{°}. However, the significant difference is that *S*(0)<0 will not preclude a physically realistic solution, because *s*^{in}∼*O*(1) will imply that it is still possible that *s*(0)>0.

Although the above constitutes one possible asymptotic structure, it is not the only one. The above assumes that there is two-phase flow in the entire GDL, but it is evident that it can break down if *θ*_{c} approaches 90^{°}. In this case, which is clearly only possible for hydrophobic GDLs if *η*>*η*^{crit}, it is natural to suppose that there would be a gas-phase only region near the catalyst region, and a two-phase region further out where *s*∼*O*(1). In this two-phase region, equation (3.9) reduces at leading order to
4.26using transformation (4.3) this time, we have
4.27where . Although no further analytical progress is possible for this case, even this much provides a useful qualitative interpretation for the forthcoming numerical work.

Lastly, note also that this solution structure also appears to be possible for the case of a hydrophilic GDL in a cell operating at low current density, i.e. for *η*<*η*^{crit}; thus, water is produced in such small quantities at the catalyst that it evaporates, but liquid water is sucked from the GDL/channel interface.

For this case, it is easier to see what is happening by considering the two-phase only case; thus, we have *Λ*=1 and, in addition, we set First, it is convenient to note that, regardless of the value of *ε*,
4.28which then gives
4.29and
4.30the boundary conditions are, at *Y* =0,
4.31
4.32and, at *Y* =1,
4.33At first sight, equation (4.29) suggests that there might be a boundary layer in *s*, because it contains a small parameter in front of the second derivative in *s*; surprisingly, perhaps, this turns out not to be the case. Instead, a consistent asymptotic structure is obtained by setting
4.34equations (4.29) and (4.30) then become, at leading order in *ε*,
4.35and
4.36where *ϱ*^{in}, , , and Δ^{in} denote the values of *ϱ*^{(g)}, *ω*^{(g)}_{O2}, and Δ, respectively, that are obtained when setting , . At leading order, boundary conditions (4.31) and (4.32) both reduce to
4.37whereas (4.33) becomes
4.38The shortfall in boundary conditions at *Y* =0 leads us to consider (4.31) and (4.32) at the next order in *ε*, which gives
4.39and
4.40this introduces the, as yet, unknown, *s*_{1} which can fortunately be eliminated to give
4.41at *Y* =0. Although this gives us the boundary condition that we need, difficulties remain because *κ*^{(g)}_{rel}=0 when *s*_{0}=1. Noting from equation (4.36) that, for all *Y* ,
4.42where *C*_{0} is a non-zero constant, we have that
4.43leading to
4.44Also, and
because *s*_{1}(0)<0,Δ^{in}<0. For example, we see that if *κ*^{(g)}_{rel}(*s*)∼(1−*s*)^{n} as then
4.45which is physically realistic only if the positive sign is taken; in this case, we have
4.46

By this stage, we have been able to decouple the equations for *s*_{0} from those for and are left with
4.47subject to
4.48however, *C*_{0} is not known, although we also have, for all *Y*,
4.49with, in particular, at *Y* =1,
4.50Were one to solve numerically for *s*_{0}, we would need to guess *C*_{0} when solving (4.47) subject to (4.48), then update the guess using (4.50); however, the more important point is that we have shown how a self-consistent boundary-value problem has been recovered in the apparently degenerate limit for *ε*≪1.

## 5. Results and discussion

Solutions to the equations given in §3*a*,*b* were computed using the finite-element software COMSOL Multiphysics; the details of the implementation follow those given by Vynnycky [24], and are therefore not repeated here. All computations were carried out using *p*^{out}=1 atm. First of all, it is necessary to determine for which model parameter values a two-phase model is even necessary. For this, we run the gas-phase only model for a range of values of *η*,*T* and RH, where RH denotes the relative humidity at the GDL/channel interface and is given by
5.1if anywhere for 0<*Y* <1 for an (*η*,*T*,RH)-combination, then the two-phase model is necessary for that combination. The results are summarized in the electronic supplementary material, figure S2, which shows the critical value of the overpotential, *η*^{crit}, as a function of cell temperature for which two-phase flow first appears for RH = 0, 30, 60, 90 per cent thus, the two-phase model is necessary for (*η*,*T*)-combinations lying above each of these curves. The fact that *η*^{crit} is an increasing function of temperature suggests that, although both and *I*_{0} both increase with temperature, it is the increase in the former that governs the appearance of these plots; otherwise, they would have been decreasing functions of temperature.

In what follows, we take
5.2as constitutive relations for the relative permeabilities, and use the expression for *p*_{c} given in equation (1.2), with *J*(*s*) given by the Udell form [31],
5.3We present first the results for the case when there is two-phase flow only, i.e. case II, followed by the cases where there is one- and two-phase flow, i.e. I and III.

### (a) Two-phase flow only: *s*^{in}=0

First, we consider the case when there is no free boundary to solve for, as this is the case corresponding to the results in Pasaogullari & Wang [8] and Song *et al.* [23]. Electronic supplementary material, figure S3 shows *s* as a function of *Y* for *θ*_{c}=0^{°},60^{°},120^{°} and 180^{°}, using equation (2.20); for this plot, *η*=0.4 V, *T*=353 K and RH = 100 per cent. Thus, as suggested earlier, this leads to the same trend as in Pasaogullari & Wang [8] and Song *et al.* [23], with saturation levels increasing as *θ*_{c} decreases from 180^{°} towards 90^{°}, but then decreasing as *θ*_{c} decreases further; whereas the trend for *θ*_{c}>90^{°} has been interpreted as indicating that increasing hydrophobicity improves the performance of the cell, because it reduces the saturation level at the catalyst layer, an explanation for the trend for *θ*_{c}<90^{°} has never been provided.

From now on, we consider only results obtained using the formulation with equation (2.19). In addition, all subsequent calculations were carried out with *η*=0.45 V and *T*=333 K, which is a combination that lies in the two-phase region for all values of RH, according to electronic supplementary material, figure S2; consequently, we do not repeat the values of *η* and *T* in every figure caption. First of all, the results for *θ*_{c}<90^{°} remain the same as those in Pasaogullari & Wang [8] and Song *et al.* [23]; however, for *θ*_{c}>90^{°}, we find instead that there is no solution at all. The physical interpretation for this is that the hydrophobic GDL does not allow any water to be present at all; note also that these results are entirely consistent with the analysis in §4*a*. Electronic supplementary material, figure S4a shows *s* as a function of *Y* for a hydrophilic GDL for four different values of *ε*. Most striking here is the limit as for which as predicted by the analysis in §4*b*. This is further reinforced by electronic supplementary material, figure S4b, which shows *x*^{(g)}_{O2} as a function of *Y* for the same values of *ε* and demonstrates that tends to its inlet value everywhere. Lastly, even though *x*^{(g)}_{O2} does not decrease from its inlet value in this limit, the current density tends to zero as as will be shown later in electronic supplementary material, figure S11, precisely because .

### (b) Two-phase flow only: *s*^{in}>0

For this case, it is possible to obtain physically realistic solutions for both hydrophobic and hydrophilic GDLs, as shown in the electronic supplementary material, figures S5 and S6, respectively, which show *s* as a function of *Y* for four different values of *ε* and *s*^{in}=0.2. From these plots, it is evident that whereas solutions can be obtained for all values of *θ*_{c} for the hydrophilic GDL, this is not the case for the hydrophobic GDL. This suggests that for any given value of *s*^{in}, there will be a critical wetting angle such that two-phase only solutions can only be obtained if *θ*_{c} is greater than this angle; we return to this point in §5*b*(ii) and in the electronic supplementary material, figure S12. Observe also that although *s*(0) will be lower for hydrophobic GDL than for hydrophilic GDL, *s*(0) decreases as *θ*_{c} decreases from 180^{°} to 90^{°}.

### (c) One- and two-phase flow: *s*^{in}=0

For this case, we need only to consider *θ*_{c}<90^{°}, because the earlier analysis in §4*a* indicates that there will be no solutions for *θ*_{c}>90^{°}. Electronic supplementary material, figure S7 shows *s* as a function of *Y* for a hydrophilic GDL for four different values of *ε* and RH = 50 per cent. The figure shows the same behaviour as electronic supplementary material, figure S4(a) for *s*(0) in the limit as (). However, the additional feature in this plot is how the location of the interface between one- and two-phase flow varies with *θ*_{c}; in general, increasing the value of *θ*_{c} shifts this towards the inlet.

Electronic supplementary material, figure S8(a) and S8(b) show *Λ* and *I*, respectively, as functions of *θ*_{c} for four values of RH. Notable here is that, apart from when *θ*_{c} approaches 90^{°}, *Λ* and *I* are relatively insensitive to *θ*_{c}. The reason for this can be seen from the governing equations. As *θ*_{c} increases from 0^{°} to 90^{°}, decreases from one to zero; however, even when *θ*_{c}≈80^{°}, i.e. only an order of magnitude decrease from its value when *θ*_{c}=0^{°}. Thus, the effect of will only really be substantial for 80^{°}≤*θ*_{c}≤90^{°}, and in this range, the effect is further magnified by the fact that multiplies the highest derivative in *s*; consequently, the profile for *s* in the electronic supplementary material, figure S7 changes dramatically as *θ*_{c} approaches 90^{°}.

### (d) One- and two-phase flow: *s*^{in}>0

For this case, we need only to consider *θ*_{c}>90^{°}, because the case for *θ*_{c}<90^{°} has already been considered in §5*b*. Electronic supplementary material, figure S9 shows *s* as a function of *Y* for a hydrophobic GDL for four different values of *ε* for *s*^{in}=0.1; in general, decreasing the value of *θ*_{c} shifts the location of the two-phase interface towards the inlet.

Electronic supplementary material, figure S10 shows *Λ* as a function of *θ*_{c} for four values of *s*^{in}. From this, it is evident that as *θ*_{c} increases, *Λ* decreases; thus, water is pushed from the inlet towards the interior of the GDL. Moreover, even comparatively small changes in *s*^{in} lead to large changes in *Λ*. Although the usual thinking is that enhancing the hydrophobicity of the GDL should cause the water to push towards the GDL/channel interface, it should be remembered that *s*^{in}>0 in this example; hence, the result is a competition between the effects of different sources of water and, in this case, the water at the GDL/channel interface appears to win out. Consequently, electronic supplementary material, figure S10 expresses the fact that, as hydrophobicity increases, water is pushed away from the dominant water source.

Electronic supplementary material, figure S11 shows *I* as a function of *θ*_{c} for the same four values of *s*^{in}; on this plot, as distinct to the electronic supplementary material, figure S10, we have also superposed results for 0^{°}<*θ*_{c}<90^{°}, in order to give a quantitative comparison of the effects of hydrophobic and hydrophilic GDLs. There are several points to note here:

— for the hydrophilic GDLs, as because but this is not the case for hydrophobic GDLs, for which ;

— the current densities obtained with hydrophobic GDLs are consistently significantly higher than for hydrophilic GDLs; and

— for hydrophobic GDLs,

*I*increases as*θ*_{c}increases, which is the only result that matches with those from the earlier hydrophobic theory, albeit for a different reason. The increase in*I*implies that*x*^{(g)}_{O2}at the catalyst layer increases as*θ*_{c}increases. Furthermore, as seen from the electronic supplementary material, figure S9, as*θ*_{c}increases, liquid water occupies more and more of the GDL; this has the effect of increasing the gas pressure gradient—although we do not present the graph here—and, hence, gas velocity in the portion of the GDL where there is one-phase flow. In turn, this leads to an enhancement in the transport of oxygen by convection, and hence helps to explain why*I*increases with*θ*_{c}. Furthermore, this also explains why increasing*s*^{in}for a fixed value of*θ*_{c}leads to an increase in*I*.

Lastly, note that we do not compute any solutions for exactly *θ*_{c}=90^{°}. Mathematically, this would not be possible because multiplies the highest derivative in *s*; instead, we are finding solutions when *θ*_{c} is just either side of 90^{°}. The reason why there is a significant difference in cell performance in the electronic supplementary material, figure S11 either side of *θ*_{c} can be traced back to the electronic supplementary material, figures S5 and S9 for *θ*_{c}<90^{°} and *θ*_{c}>90^{°}: for the former, at the catalyst layer, so that for the latter, the extent of the region of two-phase flow diminishes, but there is no reason why .

### (e) Discussion

The differences between the results of earlier work and our approach are summarized qualitatively in the electronic supplementary material, figure S12, which shows the nature of the flow as a function of the water saturation at the GDL/channel interface, *s*^{in}, and the wetting angle of the GDL, *θ*_{c}. The two approaches give identical results if *θ*_{c}<90^{°}, but vastly different results for *θ*_{c}>90^{°}:

— if

*s*^{in}=0, there is no solution at all for hydrophobic GDL;— if

*s*^{in}>0, but not too large, the catalyst layer is free of liquid water, although it is present in the GDL further towards the GDL/channel interface; and— for yet larger values of

*s*^{in}, liquid water is present everywhere in the GDL, yet the saturation level at the catalyst layer will be lower for hydrophobic GDLs than for hydrophilic GDLs.

Perhaps the starkest and most surprising result is that there is no solution at all when *s*^{in}=0 and *θ*_{c}>90^{°}; this has arisen because we have used equation (2.18), rather than equation (2.20), together with equation (1.2). The majority of the PEFC literature uses (2.20), although this flies in the face of the definition for *p*_{c} given in the more established porous-media literature. Note that replacing by in equation (1.2) and using (2.20) would lead to the same conclusion as using (1.2) and (2.18), and there may be a temptation to do this in the belief that *p*_{c} should be positive. However, we can find no physical argument why this should be the case. As counterarguments, we note the following:

— the concept of negative capillary pressure is well-known from the oil industry, and simply means that a larger water injection pressure than the oil-phase pressure has to be applied to remove oil from the sample; equation (1.2) is akin to the Young–Laplace equation, for which no absolute sign is introduced;

— Becker

*et al.*[32] obtained negative*p*_{c}via an independent route (microstructural simulations using the pore morphology method), a result that agreed well with that of Hermann & Ziegler [33]; and— Wang & Nguyen [34] obtained experimental results for

*p*_{c}which show that*p*_{c}need not be positive.

In addition, we note that Nam & Kaviany [7] use the same combination, i.e. equations (1.2) and (2.18), as we have done (see their equation 19), but then, without explanation, consider |*dp*_{c}/*ds*| instead; we can speculate that this is because they realized that there was no solution, and introduced an absolute sign to ensure that there would be. However, this has possibly led to later authors, e.g. Senn & Poulikakos [22], introducing the absolute sign, without even appealing to any belief that *p*_{c} should be positive.

## 6. Conclusions

In this paper, we have analysed in some detail a one-dimensional steady-state-generalized Darcy model for two-phase flow in the porous cathode GDL of a PEFC; the main objective was to account for an apparent anomaly in earlier modelling which implied that hydrophilic GDLs could perform as well as hydrophobic GDLs, even though experimental results do not support this conclusion. The results obtained from our re-worked theory indicate that we have indeed removed the anomaly.

Perhaps the most surprising result is that there is no solution at all when *s*^{in}=0 and *θ*_{c}>90^{°}. This comment notwithstanding, we emphasize that our results do not imply that there can never be liquid water in a hydrophobic GDL during a time-dependent process; however, they possibly indicate that water would be removed more quickly than it would in a hydrophilic GDL. This conjecture can only be resolved by a corresponding transient model, which will provide an interesting extension for this work; others are the inclusion of detailed models for: the catalyst layer and the membrane; the microporous layer; two- and three-dimensional geometries. Of particular interest will be what happens when a hydrophobic microporous layer is added to the model. For one thing, such a layer will have lower permeability because of smaller pore sizes, leading to a higher value for which has the same effect as decreasing . However, more decisive seems to be the fact that the layer is hydrophobic, because our analysis has already shown that, in steady state at least, a hydrophobic layer is less likely to retain water, leading to improved cell performance.

## Acknowledgements

The authors acknowledge the support of the Mathematics Applications Consortium for Science and Industry (www.macsi.ul.ie) funded by the Science Foundation Ireland Mathematics Initiative grant no. 06/MI/005. M.V. also acknowledges the award of a visiting researcher grant within the Government of Brazil's ‘Science without Borders’ programme.

## Footnotes

↵1 This figure, as well as all of the others for this paper, may be found online as electronic supplementary material.

- Received November 27, 2012.
- Accepted March 21, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.