F298

Journal of The Electrochemical Society, 164 (4) F298-F305 (2017)0013-4651/2017/164(4)/F298/8/$37.00 © The Electrochemical Society

Effects of Porosity and Water Saturation on the Effective
Diffusivity of a Cathode Catalyst Layer
H. Fathi,a A. Raoof,b,z S. H. Mansouri,a and M. Th. van Genuchtenb,ca Department of Mechanical Engineering, Shahid Bahonar University of Kerman, Kerman, Iranb Department of Earth Sciences, Utrecht University, Utrecht, Netherlandsc Department of Mechanical Engineering, Federal University of Rio de Janeiro, Rio de Janeiro, Brazil

The performance of a proton exchange membrane fuel cell, PEMFC, is significantly affected by the rate of oxygen diffusion throughthe cathode catalyst layer, CCL. Continuum-scale modelling of PEMFCs requires knowledge of the effective oxygen diffusivity as
a function of CCL porosity and its water saturation. To provide this functionality, we used three-dimensional pore-scale modellingto simulate the diffusion of oxygen under different liquid water saturations in CCLs having different porosity values. Solving thegoverning equations for immiscible two-phase flow, fluid distributions at different saturation levels were obtained. We show that thepresence of liquid water initiates a hindering effect by decreasing the diffusive transport of oxygen. Oxygen diffusion, including
dissolution of oxygen into the liquid water phase, was taken into account to calculate effective diffusivity values of the entire domain.
The resulting effective diffusivity values showed good agreement with values reported in the literature, which are often based on
quasi-empirical relationships. Utilizing a large number of simulation results, a correlation equation was developed for the effective
diffusivity of oxygen as a function of porosity and liquid water saturation, which is appropriate to be used for macroscopic modellingof flow and transport in CCLs.
© 2017 The Electrochemical Society. [DOI: 10.1149/2.0871704jes] All rights reserved.
Manuscript submitted September 6, 2016; revised manuscript received January 25, 2017. Published February 8, 2017.

The cathode catalyst layer, CCL, is a critical part of a proton
exchange membrane fuel cell, PEMFC, due to the loss of potential
across this section caused by transport limitations.1 CCLs have aporous structure composed of carbon, platinum, and ionomer. Aggregation of carbon particles, which are covered by platinum particles
and coated by ionomer, creates two pore categories: primary pores
within the individual agglomerates and secondary pores between the
agglomerates (Figure 1a2 ). Electrochemical reactions in the cathodeoccur at the three-phase region where ionomer, carbon, and platinum
as a catalyst for electrochemical reactions are present, together withthe gas phase. An increase in the liquid water content as a result ofelectrochemical reactions causes a decrease in the oxygen diffusivity
and lowers the active surface area.3 Predicting the cell performancerequires the coupling of the dynamics of water content with changes
in the transport properties of CCL.
A common simplifying assumption for modelling PEMFCs is toconsider the CCL as an infinitely thin layer.4–8 Therefore, the effectsof the CCL are often included as part of the imposed boundary conditions. This approach is not able to accurately consider the effect ofpotential losses, which decrease oxygen diffusion into the reaction
sites within the CCL. To obtain a more realistic model, some studies
assigned a certain thickness to the CCL.9,10 For example, to investigate combined liquid-gas flow within PEMFCs, Gurau et al.9 useda 2D model to represent the entire fuel cell, including the CCL as aseparate layer. They considered the catalyst layer to be an isotropic
porous medium and neglected the impeding effect of liquid water and,
therefore, the possible dissolution of reactants into the water phase.
Constrictions along the pathways of oxygen diffusion within the
CCL decrease the net rate of the electrochemical reactions, and hence
the current density. To investigate oxygen diffusion, some studies
considered the CCL to be a homogenous continuum.11,12 The CCL inthese studies is then assumed to contain equally-sized carbon particles
and pore spaces that are completely occupied by ionomer in which
the gas phase may dissolve. Other studies (e.g., Springer et al.1 ) assumed the catalyst layer to be a composite film of platinum/carbon
catalyst intermixed uniformly with ionomer. These various models
were able to closely match experimental data. However, experimental
observations13,14 have shown the absence of complete filling of poresby the electrolyte, which motivated several studies later to consider
agglomerated structures.15–17
Using an agglomerate model for the CCL assumes the presence ofa triple interface among pore, electrolyte, and particles of carbon supported by platinum particles. Baschuk and Li17 assumed the electrolytezE-mail: a.raoof@uu.nl

to form a thin layer around the agglomerates by covering platinum
particles as the catalyst sites. Their 2D model was later modified to
3D domains and applied to the entire cell.18 Numerical models haveapproximated the CCL by using regular19 and random20 structures. To
explore the underlying processes, pore scale simulations were applied
to the diffusion of oxygen in the CCLs of PEMFCs,21–25 however
mostly assuming very dry conditions.
Pore network modelling provides opportunity to drive macroscopictransport parameters based on pore scale processes.26,27 Raoof and
Hassanizadeh have shown that solute dispersion is a strong and nonlinear function of water saturation in granular media.28 They havefurther shown that the sensitivity of solute dispersion to change ofsaturation depends on the pore size distribution. To explore the presence of two-phase flow in the catalyst layer, some studies have usedpore-scale network modelling by assuming the existence of produced
liquid water within the variably-sized pore spaces.2,25,29 For example,
Hannach et al.2 developed a pore network consisting of spheres andcylinders representative of pore spaces to study multi-phase transportphenomena inside hydrophilic and hydrophobic CCLs. Pore networkmodels use physically-based local rules, in particular for trappingand snap-off processes, and applies simple pore geometries to increase the computational efficiency needed for modelling large domain sizes.27,30 However, other methods are required to investigate
fluid movement within pores having complex geometries. To the bestof the authors’ knowledge, there is a lack of studies at the pore scale
investigating the possible effect of produced liquid water in a CCL
on the oxygen diffusivity by applying dynamic two-phase flow simulations. A comprehensive two-phase pore-scale flow formulation isneeded to simulate the flow of immiscible liquid water under different

Figure 1. Graphical representation of (a) part of a CCL,2 (b) production ofwater around the agglomerates based on the model assumptions, and (c) thecomputational domain used for the simulations (part of the domain was removed to better visualize its internal structure).

Downloaded on 2017-02-09 to IP 131.211.208.19 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).

Journal of The Electrochemical Society, 164 (4) F298-F305 (2017)
saturation conditions, as well as for domains with different porosity
values.
In this study, we present a three-dimensional pore-scale model
to calculate the oxygen diffusivity of CCLs in the presence of liquid water. Using a finite volume method, the governing equations for
immiscible, two-phase flow were solved for domains composed ofmultiple agglomerates to obtain resulting water distributions. Simulation results were subsequently used to formulate a correlation equation for the effective oxygen diffusivity as a function of porosity andsaturation level. Finally, the accuracy of the developed relation was
validated using a number of additional independent simulations. The
developed relation is appropriate for use in continuum-scale modelling formulations where the effects of saturation and porosity on theeffective diffusivity should be provided through the use of constitutive
relations.
Model Formulation
The main objective of this study is to explore the effect of liquidwater saturation on the effective (i.e., macroscopic) oxygen diffusivity for several reconstructed CCL domains having different porosity
values. Liquid water is typically produced at the reaction sites due
to electrochemical reactions. We considered reaction sites at the surface of agglomerates by assuming a uniform distribution of platinum
particles and a high ionic conductivity of thin ionomer and carbon
all over the surface of agglomerates.2,25,29 To represent the produced
liquid water, we assigned a uniform flux of water into the secondarypores in the direction normal to the surfaces of the agglomerates. The
imposed flux causes an increase in water saturation within the CCL
pore spaces. Figure 1b shows the water distribution, schematically,
based on our model assumptions. Such a flux boundary condition implies the flow of liquid water within a hydrophilic media. Althoughhydrophobicity of a CCL is important for liquid water removal efficiency and to minimize water-covering effects at the reaction sites,
Hutzenlaub et al.31 showed that the wettability of a CCL has only alimited effect on the resulting effective diffusivity. In the model presented by Hutzenlaub et al.,31 water was assumed to fill the smallest
pore spaces prior to the larger pores for a hydrophilic domain and viceversa for hydrophobic media. They assumed that water does not exit
through domain boundaries. The flow regime, which can vary between
viscous flow or capillary flow, may be characterized using capillary
numbers and viscosity ratios of the non-wetting and wetting phases.25
The viscosity ratio in a fuel cell operating at 80◦ C is estimated to be,
M = μ2 /μ1 ∼ 18 where μ2 and μ1 are the viscosities of water and air,respectively. The capillary number, which is the ratio of viscous force
to the capillary force, can be defined as Ca = μ2 U2 /σ where U2 is thevelocity of water and σ is surface tension of water in contact with air.
The two-phase flow regime in a CCL is often controlled by capillary
forces.25,32 We have chosen a small and uniform value for the water
flux entering the domain, normal to the surface of agglomerates, tostablish such a capillary dominant flow regime. Doing so, we obtained
minimum and maximum capillary numbers of 1.1 × 10−5 and 6.7 ×
10−5 , for the smaller and the larger pores, respectively.
Flow of liquid water was obtained by direct solution of the twophase flow equations. For several water saturation levels, we performed oxygen diffusion simulations to explore the effect of water
saturation on the effective oxygen diffusivity. To obtain the distribution of water in the CCL, two equations must be solved: a momentum
equation that includes surface tension forces, and the continuity equation. The momentum equation may be written as:∂ (ρU)+ ∇. (ρU U) − ∇. (μ∇U) − ρg  



 ∂t
viscous forces

inertial forces

=−

∇p

pressure forces

−

Fs


gravity forces[1]

surfce tension forceswhere ρ is fluid density, U is the local fluid velocity, t is time, μis viscosity, p is local pressure, and Fs is the surface tension force

F299

applied at moving phase interfaces. The continuity equation is givenby:∇. U = 0[2]

When air is present in the CCL, molecular diffusion in the gas phasewill become the dominant transport mechanism since diffusion in liquid water is much slower (up to 4 orders of magnitude less) comparedto the air phase. To consider diffusion of oxygen in both the liquid andgaseous phases, the transport equation for oxygen may be written as:






∇. s DOw2 ∇cOw2+ ∇. (1 − s) DOa 2 ∇cOa 2 = 0
[3]






diffusion within the water phasediffusion within the air phase

where cOw2 and cOa 2 are oxygen concentrations in water and air, respectively, s is the local volume fraction of water, and DOw2 and DOa 2 are
the apparent diffusion coefficients [m2 .s−1 ] of oxygen in water andair, respectively.
Having two phases, Henry’s law is applied at the air-water interface
to relate concentrations of oxygen in the air and water phases:33,34
cOa 2 = HOcc2 cOw2

[4]

HOcc2

is Henry’s solubility constant. Equation 4 can be subwhere
stituted into Equation 3 to provide a governing equation in term ofcOa 2 only:s DOw2∇.
+ (1 − s) DOa 2 ∇cOa 2 = 0
[5]
HOcc2
In the absence of liquid water, all of the CCL pore space is filledwith the air containing oxygen. Under this condition, oxygen willbe transported through the medium due to the binary diffusion. Pore
spaces in the CCL are often in the range of micro- and nanoscale.
The oxygen diffusivity for this reason must account for both thebinary diffusivity of oxygen-nitrogen, DO2 −N2 , as well as the Knudsen
diffusivity, DO2 ,Kn . The dusty gas model35 may be used to obtain the
Knudsen diffusivity36 as:
DO2 ,Kn =

dp
3

8RT
πMO2

[6]

where dp is the pore diameter, R is the universal gas constant, T istemperature, and MO2 is the molecular weight of oxygen. An effective,
local scale oxygen diffusivity in air, DOa 2 , can subsequently be obtainedusing the Bosanquet equation37 as:
−1
1
1
+
[7]
DOa 2 =
DO2 −N2
DO2 ,Kn
Applying Equation 7, Knudsen diffusion is the dominant process for
relatively small pores, while bulk diffusion becomes the controlling
process for larger pores. Some studies consider oxygen diffusion to occur in both the primary and secondary pores of the CCL.21,22 However,
several studies have shown that primary pores are less accessible foroxygen and that electrochemical reactions predominantly take placeat the surface of agglomerates (i.e., at the contacts with the secondarypores), which allows one to limit diffusion only to the secondarypores.2,38,39 In this study, we hence considered oxygen diffusion to
occur only within the secondary pores. We used the average pore
diameter for calculation of the Knudsen diffusivity. Lange et al.21 calculated the local pore diameter by averaging the length of each porecomputed in thirteen different directions. They showed that similar
results for the diffusivity are obtained when considering the meanpore diameter in calculation of the Knudsen diffusivity, which is also
the choice for our study.
Numerical Simulations
Volume of fluid method.—The flow of the liquid water was simulated by applying the volume of fluid, VOF, method. Several studies
used the VOF method to solve the two-phase flow equations at the

Downloaded on 2017-02-09 to IP 131.211.208.19 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).

F300

Journal of The Electrochemical Society, 164 (4) F298-F305 (2017)

micro- and nano-scales.40–43 Bedram and Moosavi41 demonstrated theability of the VOF method to model micro- and nano-scale processes
by simulating the breakup of droplets and verifying results against
experimental observations and analytical solutions. In our study, we
solved Equations 1 and 2 in order to obtain the liquid phase volume,
Fv , in each cell of the imposed numerical mesh, such that Fv = αVcell ,where Vcell is the cell volume. The value of α will be one and zerofor completely saturated and dry cells, respectively, and between zeroand one for cells traversed by phase interfaces.44 The movement of interfaces is calculated by solving a transport equation for α as follows:∂α+ ∇. (Uα) + ∇. (α (1 − α) U r ) = 0∂t[8]in which
U r = cr |U|

∇α
|∇α|

The last term of Equation 8 is applied only to the interface. This
term represents a compression term as suggested by Weller45 in orderto minimize dispersion of the interface. The parameter cr controls
the compression of the interface. Following Weller,45 the value of cr
was taken within the range of 1 ≤ cr ≤ 4 in order to obtain a sharp
interface. The physical properties of a numerical cell were calculated
as weighted averages using the volume fraction, α, of the two fluids
occupying the cell. The fluid density at any point within the domain
is then given by:ρ = αρw + (1 − α) ρa[9]where ρw and ρa are the densities of water and air, respectively. The
viscosity, μ, is calculated in a similar way. The surface tension force,
Fs , can be obtained as:
F s = σκ (x) n[10]where n is a unit vector normal to the interface, defined as:
n=

∇α
|∇α|

[11]

and κ is the curvature of the interface, which is obtained usingκ (x) = ∇. n[12]

No-flow boundary condition was used for the flux of water across thesurrounding boundaries. It should be noted that in a real CCL, water
may exit the domain boundaries in that the amount of the remaining
water depends on pore structures and pore surface chemistry. To enable
a wide range of water saturations in this study, we applied closed
boundaries. If not done so, the maximum achievable water saturationlevel would generally be less than 0.5 (depending on the wettabilityof the domain), thus preventing a study of the effective diffusivity athigher saturation levels.
The possible impact of some of the assumptions in our model, suchas using closed boundaries, using a flux of water to represent liquid
production at the reaction surfaces, and the wettability effects, were
explored further. This was done by constructing an additional model(referred to as the ‘detailed’ model) that employs open boundaries,
water production, as well as hydrophobic surfaces. The detailed model
is described in the appendix where we show that the assumptions madein this study are appropriate and have little effect on the calculated
effective diffusivity values.
Effective oxygen diffusivity.—To obtain the effective oxygen diffusivity of the entire three-dimensional domain, oxygen diffusion was
simulated assuming different saturation levels. We invoked Dirichletboundary conditions for the oxygen concentration at the two opposite
boundaries, while the remainder of the external boundaries were considered to be symmetric. No oxygen transfer into the agglomerates
was considered. Equation 5 was solved for several partially-saturated
domains having different porosities. The distribution of liquid water in each domain, and for various saturation levels, was obtained

Table I. Parameters used in the simulations.
Property

Value

Temperature, T (K)Pressure, p (atm)
Agglomerate diameter (nm)DO2 −N2 (cm2 /s)w
DO(cm2 /s)2
HOcc2
σ (N/m)353
1
1000.27541
6.395 × 10−542
42.78543
0.06344by solving Equations 1, 2 and 8. We used the open source software,
OpenFOAM,46 to solve the governing equations. Palakurthi et al.43
verified that simulations of liquid penetration based on VOF method
using the OpenFOAM software followed Lucas–Washburn kinetics at
the micro-scale, except during the very initial stages dominated by theinertial forces.
The effective diffusivity, Deff , for each three-dimensional domain
with a given porosity and saturation was calculated as:
Deff =J. L
A (c1 − c2 )[13]

where J is the total flux of oxygen calculated over the outlet boundary
in the unit of [mol.s−1 ], L is the length of the domain, c1 and c2 are
oxygen concentrations (assigned to be one and zero at the inlet andoutlet boundaries, respectively), and A is the cross-sectional area of thebounding box of the domain perpendicular to the diffusion direction.
For each realization of the porous medium, three cases were run by c1
and c2 at opposite boundaries normal to the three orthogonal directions(x, y, and z). Approximately equal results were obtained for the threecases. However, average values of the effective diffusivities obtainedusing the three orientations are shown later along with the simulation
results.
To compare results from the various simulations with different
porosities and saturations, we defined a normalized effective diffusivity, D, as:D=

Deff
DO2 −N2[14]

in which the binary diffusivity of oxygen-nitrogen is used as the reference diffusivity. The parameters used in the simulations are shownin Table I.
Pore scale modelling.—The CCL in this study was reconstructedfollowing an algorithm proposed by Lange et al.21 After placing thefirst sphere randomly in the domain, subsequent spheres were located with a specific probability of maximum overlap with previous
spheres. Lange et al. showed that overlap values less than 30% of
the agglomerate radius had negligible effect on the effective diffusivity of the domain. We applied different levels of overlap with alimit of 20% to generate domains with different porosities. Figure 1cshows a computational domain consisting of 855 agglomerates having
a size of 100 nm and a maximum overlap of 10%. The total computational domain, shown in Figure 1c, is a cube with sides of 1 μm
and having a porosity of 0.55. A wide range of values for the diameters of the agglomerates in a CCL have been reported in numericalstudies.51 Ma et al.52 reported that carbon particles having diametersof 30–40 nm aggregate to form agglomerates of 100–120 nm in diameter. This size range is consistent with SEM images shown in an
in situ visualization of the CCL by Zhang et al.53 Similar to several
other numerical studies,54,55 we assumed agglomerates with a size of
100 nm.
To investigate the effect of finite domain size (i.e., the representative elementary volume, REV), we simulated diffusion in several
three-dimensional pore spaces by systematically increasing the domain size. Figures 2a and 2b show the effect of domain size on thesimulation results for completely dry conditions as well as for varying

Downloaded on 2017-02-09 to IP 131.211.208.19 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).

Journal of The Electrochemical Society, 164 (4) F298-F305 (2017)

F301

Table II. Simulation results for the domain consisting of
uniformly-spaced spheres, and comparison with the Bruggeman
correlation.

Figure 2. The effect of domain size on the results for different realizations
for a) dry conditions with a porosity of 0.54 and b) saturated conditions witha porosity of 0.55.

water saturations, respectively. Figure 2a shows that a domain size of
1.0 × 1.0 × 1.0 μm3 can be considered to be large, with deviations
with larger sizes being less than 1%, thus providing statistically stablediffusivity values independent of the domain size.

Number of spheres

Porosity

D

1
8
27
64
125
Bruggeman

0.47600.47400.47560.4756
0.47670.47640.33780.32930.32300.32250.32040.3288using porous domains composed of uniformly spaced, equally-sized
solid spheres.57 The porosity of such a media can be estimated as(1 − π/6), regardless of sphere diameters. This condition leads to D
= (1 – π/6)3/2 ≈ 0.3288. Table II shows calculated normalized effective diffusivities of the domain consisting of different numbers ofgrains which were uniformly spaced, as well as the normalized effective diffusivity obtained using the Bruggeman experiments. Please
notice that the simulation results are in very good agreement with
the Bruggeman experiments. Similar findings have been reported byLange et al.21 The slight deviations in the results may be partially due
to our use of a structured mesh, which cannot reproduce sphericalgeometries exactly.
We further compared our results with calculations reported byLange et al.21 They randomly placed carbon particles as equally-sized
spheres in a computational domain using the algorithm explained in
Pore scale modelling section. They obtained the effective diffusivityof the domain for different pore volume fractions and explored theeffect of Knudsen diffusion. The diameter of each carbon particle intheir simulations was 24 nm bounded by ionomer with a thickness of
4 nm in a domain with the size of 200 × 200 × 200 nm3 . For ourcomparisons we used a sphere diameter of 32 nm in a domain of 320
× 320 × 320 nm3 . The ionomer film can be treated as part of the
carbon particle because of its smaller oxygen diffusivity compared tothe binary diffusivity of oxygen-nitrogen (up to 4 orders of magnitudeless).21
Calculated normalized effective diffusivities for different porosities with and without Knudsen diffusion obtained using our model andthe model of Lange et al.21 are provided in Figure 3. Our simulation
results are in good agreement with the results from Lange et al., withor without including the Knudsen effect. Figure 3 also reveals the

Results and Discussion
In this section, we provide first results of the effective diffusivity
for dry conditions since this situation has been explored in several previous studies.21,39,56 Subsequently, the effect of water saturation on theeffective diffusivity of CCLs having different porosities is explored
for the purpose of formulating a possible correlation equation. Next,
the importance of considering immiscible two-phase flow for simulating a CCL is shown by comparing our results with those obtained
assuming miscible flow of water into the gas phase. At the end, theeffect of agglomerate sizes on the effective diffusivity of the CCL isinvestigated.
Diffusivity under dry conditions.—Simulation results of the diffusivity are often validated against Bruggeman correlation. Bruggeman
relation57 describes tortuosity, τ, as τ = ε1/2 , where ε is the porosity.
Tortuosity may be correlated to the normalized effective diffusivity, D,
through the expression58 D = ε/τ. Substituting the Bruggeman relation
into the effective diffusivity expression provides D = ε3/2 . This classical equation is known as Bruggeman correlation for the effective
diffusivity.21 Moreover, Bruggeman conducted several experiments

Figure 3. Comparison of normalized effective diffusivity values obtained in
the present study with results reported by Lange et al.21 with and without
Knudsen diffusion, together with the Bruggeman correlation.57

Downloaded on 2017-02-09 to IP 131.211.208.19 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).

F302

Journal of The Electrochemical Society, 164 (4) F298-F305 (2017)

Figure 4. Liquid water distributions in the domain at saturation values of (a)0.2, (b) 0.4, (c) 0.6, and (d) 0.8. The porosity of the domain was 0.61. Liquidwater is shown in blue and the gas phase in white.

considerable effect of Knudsen diffusion which is due to the presenceof small-sized pores.
Diffusivity under different saturation levels and porosity
values.—We used the three-dimensional water distributions resulting from the two-phase flow calculations at different saturation levels
to simulate the diffusion of oxygen. Note that Figure 2 showed thatthe chosen domain size of 1.0 × 1.0 × 1.0 μm3 was large enoughfor simulations under dry and saturated conditions. Figure 4 shows
the liquid water distributions at several saturation levels in a domain
having a porosity of 0.61. The plots indicate that liquid water invades the pore spaces between agglomerates due to capillary forces
to change local water saturations within the individual pores. Figure
5b shows the steady state distribution of oxygen within the domain
having a porosity of 0.61 and saturation level of 0.2 (shown in Figure5a). Clearly, for the applied oxygen concentrations at the inlet andoutlet boundaries, higher water saturations lead to lower local rates ofoxygen diffusion due to the significantly lower diffusion coefficient ofoxygen in liquid water compared to its value in the air phase. Figure 6
further shows how for several different porosities the stablished fluxof oxygen decreases as a function of liquid saturation.
Analyzing the results from all simulations (i.e., for different porosities and saturation values), we developed the following correlation toestimate the normalized effective diffusivity, D, as a function of porosity, ε, and the average water saturation, S, of the entire domain:


S
3.2
[15]
D = ε(2.6+4ε ) (1 − S)(0.54+1.4 ) + βS

Figure 5. Distributions of the simulated (a) liquid water distribution (withliquid water shown in blue and the gas phase in white) and (b) oxygen concentration distribution within the computational domain having a porosity valueof 0.61 and a saturation level of 0.2.

of the fitted equation, we used Equation 15 to predict effective diffusivities at porosities of 0.46 and 0.68 which are, respectively, insideand outside of the parameter range used to develop this equation. We
compared the results obtained using the pore scale model for two
domains having different porosities with the values obtained using
Equation 15. Figure 8 verifies the accuracy of Equation 15 to provide
effective diffusivity estimates for a wide range of porosity values.
Given the significantly smaller diffusion coefficient of oxygen inwater compared to the air phase, the effect of increasing liquid wateris very similar to the effect of a decreasing porosity.59 This implies
that the effective diffusivity may be reduced by a factor (1-S)1.5 , similarly as the effect of porosity in the original Bruggeman correlation.
To evaluate this assumption, Figure 9 shows the dependency of normalized gas diffusion coefficient on water saturation using our results
as well as using a Bruggeman type correction. Figure 9 indicates thatusing the Bruggeman type correction provides a good match of theresults of the pore scale simulations on immiscible flow. The Bruggeman type correction hence could be used to account for the effectof water saturation on the effective diffusivity for immiscible flow

where
β=

DOw2
DO2 −N2.
3.2

The second term on the right-hand side of Equation 15, i.e.,βε(2.6+4ε ) ,is for oxygen diffusion in a fully-saturated domain. This term provides negligible contribution when liquid water starts to block diffusion through the gaseous pathway in the domain. Figure 7 provides
a three-dimensional surface representation of Equation 15. This figure suggests a monotonic increase in the effective diffusivity withan increase in air content (i.e., with decreasing water saturation) andporosity. Equation 15 was developed using a total of 44 simulations
with porosity values between 0.34 and 0.61. To validate the accuracy

Figure 6. Relationship between the calculated flux of oxygen per unit areaand liquid water saturation for different porosities.

Downloaded on 2017-02-09 to IP 131.211.208.19 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).

Journal of The Electrochemical Society, 164 (4) F298-F305 (2017)

Figure 7. Surface graph of Equation 15 showing variations in the normalized
effective diffusivity as a function of liquid saturation and porosity.

conditions, except perhaps for a slight overestimation of the effective
diffusivity at relatively high saturation levels.
Immiscible two-phase flow effects.—In this study, the flow of liquid water was simulated by solving the governing two-phase flow
equations to capture as best as possible the immiscible distributionof water and air phases, including their well-defined interfaces due tocapillary effects. However, in order to lower computational times andnumerical complexities, one may assume the miscible displacementof water in the gas phase (i.e., single-phase flow conditions in contrast to the two-phase flow conditions). This implies that for a given
saturation, the water will be uniformly distributed within the domain.
To evaluate this assumption, we compared the two methods (i.e., simulations considering miscible and immiscible water displacement).
Figure 9 shows the discrepancy between the two methods in termsof calculated effective diffusivity values. Results indicate that the assumption of miscible flow significantly overestimates the normalized
effective diffusivity. This is due to the fact that during two-phase flow,
capillary forces cause water to invade into the pore spaces between
the agglomerates. Liquid water then effectively blocks the pore spaceconnectivity towards the air phase, which contributes the most to
oxygen diffusion. Significantly lower effective diffusivity values arecaused due to the absence of a percolating air pathway through thedomain. Therefore, the capillarity of the liquid water phase is a key
process which should be accounted for by solving the two-phase flow
formulations.

Figure 8. Validation of the accuracy of the correlation given by Equation 15
for porosity values, ε, of 0.46 and 0.68 within and outside of the range ofporosities used to derive Equation 15, respectively.

F303

Figure 9. Effect of water saturation on the normalized effective diffusivity as
obtained by using the Bruggeman type correction59 as well as the assumptionsof immiscible and miscible flow in the porous medium. Results are for a domain
having a porosity of 0.55.

Diffusivity for different agglomerate sizes.—In the adopted porescale modelling approach, the size of the solid agglomerates determines the size of the resulting pore spaces which, in turn, control the
Knudsen diffusivity through Equation 6. Figure 10 shows the pore size
distributions for three domains having the same total porosity valueof 0.55, but composed of spheres with three different diameters, toproduce differently-sized pore spaces. An increase in the agglomerate size apparently creates a much wider pore size distribution, whichleads to a higher mean pore diameter and hence lower Knudsen effect.
In contrast, a decrease in agglomerate size decreases the mean pore
diameter, leading to a higher Knudsen effect. For the same domains,
Figure 11 shows the effect of agglomerate size on the normalized effective diffusivity. The domain with the largest spheres produced thehighest values of the normalized effective diffusivity due to decreasing
Knudsen effects.
Conclusions
In this study we used a three-dimensional pore scale model toinvestigate diffusive transport of oxygen through a cathode catalyst
layer, CCL. A large number of simulations were conducted using domains having different porosity values and a wide range of liquid water
saturations. Water distributions within the CCL were obtained using
the volume of fluid method for solving the governing two-phase flow
equations. Effective diffusivity values were calculated based on the

Figure 10. The effect of agglomerate size on the pore size distribution ofa domain with the porosity of 0.55. The plot shows the probability density
function (PDF) using Weibull distibution.60 The coefficient of variation of
three pore size distributions are equal.

Downloaded on 2017-02-09 to IP 131.211.208.19 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).

F304

Journal of The Electrochemical Society, 164 (4) F298-F305 (2017)

Figure 11. Effect of agglomerate size on the normalized effective diffusivity
for different saturation levels. The size of the cubic domain for each realization
was ten times the agglomerate diameter, while the porosity was 0.55. Data
points were connected to improve visibility.

diffusion of oxygen through the resulting two-phase air-water system
considering oxygen dissolution in the water phase as well as Knudsen
effects. The resulting effective diffusivity values for different porosity
and saturation levels were compared with the Bruggeman empiricalcorrelation, as well as with numerical results reported by Lange etal.21 We found that choosing 100 nm size agglomerates produced asystem with significant effect of Knudsen diffusion. Based on a large
number of simulations, a correlation equation was developed to obtain the effective diffusivity of oxygen as a function of porosity andliquid water saturation. This constitutive relation appears promising
for macroscopic continuum scale numerical studies demanding valuesof the oxygen diffusivity as a function of the porosity and saturation
state of a CCL.
Our results indicate a significant decrease in the total oxygen flux
across the domain at high saturation levels. This reduction is due tothe impeding effects of the produced liquid water within the porous
medium. This produced water tends to invade the pore spaces within
the CCL, leading to fewer percolating and more tortuous pathways foroxygen diffusion through the domain. A decrease in porosity resulted
in very similar effects on the oxygen diffusivity as an increase in theliquid water saturation.
We showed that the assumption of miscible displacement of water will significantly overestimate calculated diffusivity values. In ourstudy we used an equivalent flux at the surface of the agglomerates to account for water production by electrochemical reactions.
Future studies are required which directly implement the electrochemical reactions in the two-phase flow simulations, including thepossibility of ionomer coatings causing wettability changes. In thisstudy we neglected any evaporation of liquid water within the CCL.
This assumption is expected to be appropriate since the gas phase in
the pore spaces becomes saturated with water vapor in a short time,
which causes the net evaporation rate of water to become zero suchthat water transport is dominated by capillary flow.61 Although ourstudy applied three-dimensional simulations, the use of real CCL pore
structures (i.e. those obtained using non-destructive imaging techniques), should provide a more realistic representation of the solid
surfaces.

Figure A1. A schematic of the numerical cells used in the computational
procedure.

Due to the presence of micro- and nanoscale pores in a CCL, and the strongly coupled two-phase flow and water production processes involved, solving the two-phase flowat each time step of water production can be extremely time consuming. Additionally,
numerical stability requires very small time steps. For this reason we used a sequential
approach which assumes that cells adjacent to the surface of agglomerates are saturated(i.e., the water fraction, α, is set to a value of 1), after which the two-phase flow equations are solved to lead to the equilibrium water distribution. During the two-phase flow
simulation, water will move within the domain, resulting in desaturation of some cells inthe domain, in particular those cells which are adjacent to the hydrophobic surfaces of the

Appendix: Description of the Detailed Model
This section describes the more detailed model we constructed to explore the validityof assumptions made in this study. The detailed model employs open boundaries, waterproduction at the reaction surfaces as well as hydrophobic surfaces. The model assumes
that when the molar ratio of liquid water in a grid cell of the computational domainreaches a value of one, the local numerical cell is considered to be water saturated.
Figure A1 provides a schematic part of the numerical domain based on the detailed model
assumptions. The figure shows that cells adjacent to the surface of agglomerates form thereaction zone where liquid water is being produced.

Figure A2. Simulated water distributions for a hydrophobic domain with acontact angle of 150 and a porosity of 0.55. The figure shows the domain
having an initial saturation of 0.20 (a) and a final saturation of 0.28 (b). Liquidwater is shown in blue.

Downloaded on 2017-02-09 to IP 131.211.208.19 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).

Journal of The Electrochemical Society, 164 (4) F298-F305 (2017)

Figure A3. Effective diffusivity versus water saturation as obtained by thepresent (standard) model with closed boundaries and the detailed model whichallows water to exit the domain. The final water saturation of the domain with
open boundaries was 0.28 when the contact angle was 150.

agglomerates. This procedure of placing water saturated cells next to the agglomerate surfaces and solving the two-phase flow equations is repeated until reaching an equilibrium
water distribution within the whole domain. For the boundary conditions in the two-phaseflow simulations we assumed the surface of the agglomerates to be solid hydrophobic
surfaces with a contact angle of 150, while pressure outlet conditions were applied to theremaining boundaries. When the water flow system reached steady-state conditions, thefinal saturation value was found to be close to 0.28 for the domain having a contact angleof 150 and a porosity of 0.55 (Figure A2).
We compared the normalized effective diffusivity values of the two models (i.e., withand without letting water exit the domain) for the same saturations. The final saturation(at which steady state flow was reached) for the CCL with water removal was 0.28 forthe domain with contact angle of 150. A comparison with results for the case with closed
boundaries is shown in Figure A3. Although there is some difference between the twocases, the general agreement is good and hence justifies the use of closed boundaries toenable our simulations to cover a wide range of water saturations.
The detailed model represents more physics including: the hydrophobic CCL surfaces
as well as water production at the reaction sites instead of assigning a flux of water. The
relatively good agreement between results obtained with the standard and detailed models
shows that the assumption of no water exiting the domain, using the water flux fromthe surface of the agglomerates, and assuming a hydrophilic domain, provides valuable
insight into the effective diffusivity of CCLs. We note that the simulation time for thedetailed model was significantly longer (approximately 10 times), which made it very
difficult for us to run a large number of simulations as needed for this study.

References
1. T. E. Springer, M. S. Wilson, and S. Gottesfeld, J. Electrochem. Soc., 140 (12), 3513
(1993).
2. M. El Hannach, J. Pauchet, and M. Prat, Electrochim. Acta, 56 (28), 10796 (2011).
3. F. Barbir, PEM fuel cells: Theory and practice, Academic Press, Storrs, United States(2012).
4. D. M. Bernardi, J. Electrochem. Soc., 137, 3344 (1990).
5. T. E. Springer, T. A. Zawodzinsky, and S. Gottesfeld, J. Electrochem. Soc., 138 2334
(1991).
6. T. V. Nguyen and R. E. White, J. Electrochem. Soc., 140 (8), 2178 (1993).
7. J. H. Nam and M. Kaviany, Int. J. Heat Mass Transfer, 46 (24), 4595 (2003).
8. T. Berning and N. Djilali, J. Electrochem. Soc., 150 (12), A1589 (2003).
9. V. Gurau, H. Liu, and S. Kakac, AlChE J., 44 (11), 2410 (1998).
10. M. Tohidi, S. Mansouri, and H. Amiri, Int. J. Hydrogen Energy, 35 (17), 9338 (2010).

F305

11. D. M. Bernardi and M. W. Verbrugge, AlChE J., 37 (8), 1151 (1991).
12. D. M. Bernardi and M. W. Verbrugge, J. Electrochem. Soc., 139 (9), 2477
(1992).
13. E. Ticianelli, J. Beery, and S. Srinivasan, J. Appl. Electrochem., 21 (7), 597(1991).
14. M. Uchida, Y. Aoyama, N. Eda, and A. Ohta, J. Electrochem. Soc., 142 (12), 4143
(1995).
15. K. Broka and P. Ekdunge, J. Appl. Electrochem., 27 (3), 281 (1997).
16. M. Eikerling and A. Kornyshev, J. Electroanal. Chem., 453 (1), 89 (1998).
17. J. Baschuk and X. Li, J. Power Sources, 86 (1), 181 (2000).
18. D. H. Schwarz and N. Djilali, J. Electrochem. Soc., 154 (11), B1167 (2007).
19. G. Wang, P. P. Mukherjee, and C.-Y. Wang, Electrochim. Acta, 51 (15), 3139
(2006).
20. G. Wang, P. Mukherjee, and C. Wang, Electrochim. Acta, 51 (15), 3151 (2006).
21. K. J. Lange, P.-C. Sui, and N. Djilali, J. Electrochem. Soc., 157 (10), B1434(2010).
22. K. J. Lange, P.-C. Sui, and N. Djilali, J. Power Sources, 208 354 (2012).
23. R. Singh, A. Akhgar, P. Sui, K. Lange, and N. Djilali, J. Electrochem. Soc., 161 (4),F415 (2014).
24. P. P. Mukherjee and C.-Y. Wang, ECS Trans., 3 (1), 1085 (2006).
25. P. P. Mukherjee, Q. Kang, and C.-Y. Wang, Energy Environ. Sci., 4 (2), 346(2011).
26. L. Vasilyev, A. Raoof, and J. M. Nordbotten, Transport Porous Med., 1 (2012).
27. A. Raoof, H. Nick, S. Hassanizadeh, and C. Spiers, Comput. Geosci., 61 160(2013).
28. A. Raoof and S. Hassanizadeh, Water Resour. Res., 49 (4), 1943 (2013).
29. R. Alink and D. Gerteisen, Energies, 6 (9), 4508 (2013).
30. A. Raoof and S. M. Hassanizadeh, Transport Porous Med., 81 (3), 391 (2010).
31. T. Hutzenlaub, J. Becker, R. Zengerle, and S. Thiele, J. Power Sources, 227 260(2013).
32. P. P. Mukherjee, C.-Y. Wang, and Q. Kang, Electrochim. Acta, 54 (27), 6861
(2009).
33. M. G. Mahmoodlu, N. Hartog, S. M. Hassanizadeh, and A. Raoof, Chemosphere, 91(11), 1534 (2013).
34. M. G. Mahmoodlu, S. M. Hassanizadeh, N. Hartog, and A. Raoof, J. Contam. Hydrol.,164 193 (2014).
35. E. A. Mason and A. Malinauskas, Gas transport in porous media: the dusty-gas
model, Elsevier Science Ltd (1983).
36. M. M. Mench, Fuel Cell Engines, 285 (2008).
37. W. Pollard and R. D. Present, Phys. Rev., 73 (7), 762 (1948).
38. R. Wu, Q. Liao, X. Zhu, and H. Wang, Int. J. Hydrogen Energy, 37 (15), 11255(2012).
39. S. Litster, W. Epting, E. Wargo, S. Kalidindi, and E. Kumbur, Fuel Cells, 13 (5), 935(2013).
40. D. Banerjee, NSTI Nanotech., 10 (2005).
41. A. Bedram and A. Moosavi, J. Appl. Fluid Mech., 6 (1), 81 (2013).
42. A. Bedram, A. Moosavi, and S. K. Hannani, Eur. Phys. J. E, 38 (9), 1 (2015).
43. N. K. Palakurthi, S. Konangi, U. Ghia, and K. Comer, Int. J. Multiph. Flow, 77 48(2015).
44. C. W. Hirt and B. D. Nichols, J. Comput. Phys., 39 (1), 201 (1981).
45. H. Weller, OpenCFD Ltd., Report TR/HGW/04, (2008).
46. The OpenFOAM Foundation. http://www.openfoam.org/, (accessed 14.08.16).
47. E. N. Fuller, P. D. Schettler, and J. C. Giddings, Ind. Eng. Chem., 58 (5), 18 (1966).
48. W. Hayduk and H. Laudie, AlChE J., 20 (3), 611 (1974).
49. J. C. Amphlett, R. Baumert, R. F. Mann, B. A. Peppley, P. R. Roberge, andT. J. Harris, J. Electrochem. Soc., 142 (1), 1 (1995).
50. N. Vargaftik, B. Volkov, and L. Voljak, J. Phys. Chem. Ref. Data, 12 (3), 817(1983).
51. P. Dobson, C. Lei, T. Navessin, and M. Secanell, J. Electrochem. Soc., 159 (5), B514
(2012).
52. S. Ma, C.-H. Solterbeck, M. Odgaard, and E. Skou, Appl. Phys. A, 96 (3), 581 (2009).
53. F.-Y. Zhang, D. Spernjak, A. K. Prasad, and S. G. Advani, J. Electrochem. Soc., 154(11), B1152 (2007).
54. Q. Wang, M. Eikerling, D. Song, and Z. Liu, J. Electroanal. Chem., 573 (1), 61(2004).
55. T. Sousa and C. Rangel, J. Solid State Electrochem., 20 (2), 541 (2016).
56. K. J. Lange, P.-C. Sui, and N. Djilali, J. Power Sources, 196 (6), 3195 (2011).
57. V. D. Bruggeman, Ann. Phys., 416 (7), 636 (1935).
58. E. E. Petersen, AIChE J., 4 (3), 343 (1958).
59. A. Z. Weber and J. Newman, Chem. Rev., 104 (10), 4679 (2004).
60. W. Weibull, J. Appl. Mech., 103 293 (1951).
61. G. Lin, W. He, and T. Van Nguyen, J. Electrochem. Soc., 151 (12), A1999 (2004).

Downloaded on 2017-02-09 to IP 131.211.208.19 address. Redistribution subject to ECS terms of use (see ecsdl.org/site/terms_use) unless CC License in place (see abstract).