An Improved Agglomerate Model for the PEM Catalyst Layer with Accurate Effective Surface Area Calculation Based on the Sphere-Packing Approach

The complex porous structure of the PEM fuel cell catalyst layer (CL) necessitates the use of multiscale modeling strategies such as the agglomerate approach. In this study a 2D steady-state model for the cathode CL is developed using the spherical agglomerate approach. A new, more accurate, method is introduced to determine the effective agglomerate surface area, which plays a key role in estimating diffusion losses in the CL. Specifically, the reduction in the effective surface area due to overlapping particles is modeled geometrically based on a sphere-packing approach. In addition, the equations for the agglomerate model are reformulated to correctly account for the agglomerate surface area reduction due to overlapping particles. The importance of an accurate geometric model for the effective surface area is demonstrated by investigating the effect of CL composition on performance, and the results show that the new method provides more realistic predictions than the existing approach. Results from the new approach for optimal Pt loading, ionomer loading, and Pt |C ratio show good agreement with experimental results. © The Author(s) 2014. Published by ECS. This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 License (CC BY, http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse of the work in any medium, provided the original work is properly cited. [DOI: 10.1149/2.116406jes] All rights reserved.

Great strides have been made in proton exchange membrane (PEM) fuel cell research since the introduction of the technology in the early 1960s. Computational models represent a major tool for the design and development of PEM fuel cells because they decrease the dependence on expensive and time-consuming experimental procedures. Accurate computational models can also provide key insights into the fundamental processes that govern fuel cell performance. The overall accuracy of computational models depends strongly on the description of the catalyst layer (CL) as all the critical electrochemical reactions, and heat and mass transport processes occur within the CL. Therefore, developing better modeling strategies for the CL has been an important research objective to obtain a realistic representation of PEM fuel cell behavior.
Catalyst layer modeling approaches, which have grown in complexity in the past decades, were described in detail in. 1 Here, we present only a brief review. Among the models reported, the interface model [2][3][4] is the simplest approach as it represents all of the CL activity as boundary conditions applied at the interface between the membrane and the gas diffusion layer (GDL). The macro-homogenous model [4][5][6][7] is a more advanced approach that utilizes the effective medium theory and approximates the CL as a porous matrix of gas pores, platinum, carbon, and electrolyte. The most comprehensive device-level model to date is the agglomerate approach [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22] as it represents CL activity as a multiscale problem. Overall, the CL is assumed to consist of gas pores and aggregates of C|Pt particles that are covered by an electrolyte film. The agglomerate approach attempts to model the C|Pt aggregates (shown schematically in Fig. 1a) that are typically seen in microscopic images of the CL. For simplicity, the aggregate microstructure, denoted as an agglomerate, is idealized as in Fig. 1b to resolve the complex CL activity. Here, the agglomerate core is assumed to consist of a homogenous mixture of Pt, C, and ionomer, and is surrounded by an ionomer film of specified thickness. Essentially, the mass transport and reaction processes are solved within the agglomerate under a constant overpotential assumption, and the results are subsequently integrated with the macroscale simulation. Because it is possible to solve for the mass transport within the idealized agglomerate analytically, the proposed microstructure can be conveniently implemented in such a multiscale approach.
Recently, some researchers have proposed modifications to the classical agglomerate model by incorporating numerical solutions of the microscale problem within device-level simulations. For example, z E-mail: prasad@udel.edu Kamarajugadda and Mazumder 21 proposed a generalized geometry for the agglomerate by considering the agglomerate as two overlapped spheres with unequal radii. They incorporated the numerical solution for the overlapped agglomerates as a sub-grid-scale model into their 2D PEM fuel cell model by considering the effect of porosity. They concluded that the shape of the agglomerate only matters in the case of large agglomerates. 21 Cetinbas et al. 22 proposed a modified agglomerate approach in which the agglomerate core contained a random distribution of discrete catalyst particles instead of the idealized homogenous agglomerate mixture assumption. They showed that their modified approach was superior to the classical homogenous mixture assumption in predicting the effect of Pt loading and other composition variations within the agglomerate core. Next, they coupled their modified particle approach with a 3D cathode model 1 based on the effective surface area of the agglomerates, and obtained device-level results for the 3D PEM fuel cell cathode. In addition to these agglomerate model studies, inexpensive computational power has made it feasible to study the actual microstructure of the CL. [23][24][25][26][27][28] By incorporating the complex morphology of the CL, such models can help to characterize the CL, for example, by determining effective transport properties, which is quite difficult with experiments.
A key requirement for CL models is that they should be able to capture the effects of CL composition variations. In, 9,13,14 it was recognized that the agglomerate model is unable to realistically predict the effect of catalyst loading. All of these studies reported identical limiting currents despite wide variations in the catalyst loading. On the other hand, experimental results [29][30][31] indicate that Pt loading strongly influences fuel cell performance and the limiting current density. This discrepancy was successfully addressed in 22 by replacing the homogenous mixture within the agglomerate core with discrete catalyst particles. Consequently, it was reported that the classical homogenous assumption results in a highly reactive agglomerate core wherein the diffusion resistance imposed by the ionomer film determines the overall rate of the reaction.
Instead of modifying the agglomerate microstructure, 12,15 proposed relations to link the microstructure to macroscale volume fractions of the various CL components to improve the predictions of the classical agglomerate model. Yin 12 employed relationships which caused the ionomer film thickness δ to decrease as Pt|C was increased while maintaining a constant CL porosity v , agglomerate radius r agg , and ionomer volume fraction ε agg . In Yin's study, 12 δ can approach zero such that the agglomerate model behaves like the macro-homogenous model, so it fails to predict the limiting current in such cases. Moreover, contrary to Yin's 12 assumption, it is reported that the ionomer  fills the pores between agglomerates (the secondary pores) first, 32 implying that δ = 0. Secanell et al. 15 incorporated changes in the CL structure by varying the number of agglomerates per unit volume which concurrently changes the total specific surface area of the agglomerates with respect to the CL volume. We will denote this quantity henceforth as a agg , the effective surface area. Secanell et al. used constant agglomerate parameters and a constant CL thickness while varying the number of agglomerates to relate the CL volume fractions with the agglomerate parameters. They then proposed an expression for the effective surface area as a function of porosity and the number of agglomerates per unit volume. Here it is important to note that they modeled the reduction of the effective surface area as an arbitrary function of the CL porosity by only considering two limiting cases of CL porosity. Their relations were employed in 16,20 to study the effect of CL composition on performance.
The main goal of this study is to build a physical model to analytically quantify the reduction in the effective surface area (a agg ) as the CL porosity decreases. We have adapted the sphere-packing approach used by Taylor et al. 33 who studied the effect of bio-film growth over rigid particles on porosity and specific surface area. We developed a new method to obtain the effective surface area for two different sphere-packing schemes and the CL porosity which is calculated using model inputs. Unlike, 15 ionomer weight fraction is used as an input, and the ionomer volume fraction inside the agglomerate is calculated based on the number of agglomerates per unit volume. Furthermore, we have reformulated the agglomerate model equations to correctly account for the reduction of the effective surface area in the case of overlapping agglomerates. Next, we implemented the agglomerate approach in a 2D cathode domain to simulate the effect of ionomer weight fraction and Pt loading on performance. The influence of our effective surface area calculation method on the effect of Pt loading, ionomer weight fraction, and Pt|C ratio on performance was studied for both sphere-packing schemes and compared with the method in. 15 In addition, the predictions of our model for various CL compositions were compared against experimental data in the literature.

Model Description
Cathode model.-The 2D, steady-state cathode model accounts for the transport of chemical species, electrons, and ions in the GDL and the CL. The spherical agglomerate approach is used as an embedded model to describe the CL activity. For simplicity, the model assumes an isothermal, isobaric, and steady-state process. The computational domain encompasses half of the channel and half of the land area as illustrated in Fig. 2 and the associated geometric parameters are listed in Table I. The diffusion of oxygen, water vapor, and nitrogen in the cathode is described by the Maxwell-Stefan equations after neglecting the convective terms.
Eq. 1 is solved for the mass fraction (ω j ) of each chemical species where R j is the source term calculated from the electrochemical reactions, and mixture density ρ is calculated using the ideal gas law.
Charge transport by the electron and ion phases is governed by: where σ ef f m and σ ef f s are the effective ionic and electronic conductivities, ∅ m and ∅ s are ionic and electronic potentials, and I C L v is the volumetric current density. For both charge and mass transport the effective properties are calculated using the Bruggeman approximation.
where P is the property and i is the volume fraction of the i th phase in the porous medium. The local ionic conductivity (σ m ) is determined by using the empirical relation reported by Springer et al. 34 for Nafion.
where T is the operating temperature in K, and λ is the local water content which is expressed as a function of the local relative humidity φ R H as: 13 The source terms in Eq. 1 are only non-zero in the CL domain. The O 2 and H 2 O source terms (R O 2 and R H 2 O ) are calculated from the volumetric current generation using Faraday's law: where α O D is the electroosmotic drag coefficient which can be expressed as: 13 [9] The nominal cathode overpotential (η NC O ) is defined as the measure of total losses (including activation and ohmic losses) in the cathode by: 13 [10] where the electron potential at the current collector and the ion potential at the membrane-CL interface are the boundary conditions for the charge conservation equations. The spherical agglomerate model is used to describe the relevant electrochemical phenomena in the CL. At the microscale, the reaction rate (which is directly proportional to the current generation) can be calculated in two ways: (1) using the reactant flux at the agglomerate surface, or (2) using the effectiveness factor approach. 35 Therefore, the following equation can be written for the current generation in a single agglomerate: where r agg is the agglomerate radius, δ is the ionomer film thickness, A agg is the specific surface area of an agglomerate (ratio of agglomerate outer surface area to its volume), D is the diffusivity of oxygen in ionomer C 0 , is the oxygen concentration on the ionomer film surface, C s is the oxygen concentration on the agglomerate core surface, k c is the reaction rate coefficient, and E r is the effectiveness factor. The effectiveness factor is defined as the ratio of the actual reaction rate to the expected reaction rate if all the catalyst particles inside the agglomerate are subjected to the surface concentration of the reactant at the wall of the agglomerate. 35 Knowing the current generation for a single agglomerate, the volumetric current for the entire CL can be easily obtained. For the flux approach, the effective surface area (a agg ) must be defined as the ratio of the total surface area of all agglomerates to the CL volume. Then, by multiplying the left hand side of Eq. 11 by a agg /A agg , the current generation per unit volume of the CL can be written in terms of the flux at the agglomerate surface as in left hand side of Eq. 12. Similarly, by considering the ratio of the total volume of all agglomerates in the CL to that of a single agglomerate, the current generation per unit CL volume can be written as in the right hand side of Eq. 12.
where v is the catalyst layer porosity. Here, it should be noted that without the correction term ( ragg ragg +δ ) 3 , the right hand side of Eq. 11 gives the reaction rate per unit volume of agglomerate core rather than per unit volume of the entire agglomerate (core plus ionomer film). This correction was not taken into account in. 13,15,20 More crucially, since ( only for non-overlapping agglomerates, Eq. 12 is only valid for the case of non-overlapping agglomerates. For overlapping agglomerates, the effectiveness factor must be derived for each overlapping geometry, i.e. a agg must be computed accurately. The major contribution of this paper is to correctly estimate a agg using a sphere-packing approach, and highlight differences with previous approaches such as. 15 For the overlapping case, the effectiveness factor derived for a single agglomerate is no longer valid. In contrast, the flux approach is more suitable because it can account for the surface area reduction due to overlapping agglomerates. Assuming that C s does not change between a single agglomerate and overlapping ones, Eqs. 13-15 can be used to calculate the volumetric current generation in the CL: T h is the dimensionless Thiele modulus defined as the ratio of the surface reaction rate to the diffusion rate through the catalyst pellet: 35 D ef f is the effective diffusivity of oxygen inside the agglomerate and is calculated using Bruggeman approximation. The reaction rate coefficient depends on the fuel cell reaction kinetics according to: 15 where i re f 0 is the reference exchange current density, C O 2 ,re f is the reference oxygen concentration, α c is the cathode charge transfer coefficient, and a Pt is the total reaction area per unit volume of the agglomerate which is expressed as: where m Pt is the catalyst loading (kg/m 2 ), t C L is the catalyst layer thickness, and A 0 is the catalyst surface area per unit mass. The following correlation for A 0 was developed by 15 as a function of Pt|C ratio based on the experimental data of: 39 A 0 = 2.2779 × 10 6 (Pt|C) 3 − 1.5857 × 10 6 (Pt|C) 2 − 2.0153 × 10 6 (Pt|C) + 1.5950 × 10 6 cm 2 /g [19] Table II lists all the parameter values employed in the current simulations. The literature indicates two Tafel slopes corresponding to high and low current regimes. [36][37][38] Accordingly, we used two distinct values for i re f 0 and α c as listed in Table II. In addition to the various material property values and operating conditions, Table II also gives the agglomerate dimensions selected for this study.
All agglomerate models require the use of tuning parameters to match experimental data; in our case, we tuned the exchange current density and the agglomerate dimensions to fit the experimental polarization data of 31 as shown subsequently in Section 3. The literature reports a wide range of values for the agglomerate dimensions. For example, one of the pioneering agglomerate modeling studies employed an agglomerate radius ranging from 1-5 μm based on their SEM images. 8 On the other hand, a recent experimental study 40 suggests that the carbon particles form agglomerates of about 100-120 nm. In particular, modeling studies employ a wider range of agglomerate radius of between 0.5-5 μm compared to experimental observations such as. 40 This discrepancy can be understood by noting that the agglomerate model does not attempt to represent the actual physical microstructure of the CL. Rather, the agglomerate dimensions employed in models are selected to represent the mean path that dissolved oxygen must travel before reaching the reaction sites and thereby accurately model the various transport losses in the CL. 41 Consequently, the agglomerate sizes used in models are larger than those observed in SEM or TEM images. In our study, we set r agg = 3μm and δ = 145μm as listed in Table II.
Neglecting the losses in the anode and the membrane, the polarization curves are plotted as a function of the cathode potential which is given as: 13 Referring to Fig. 2, the boundary conditions employed to solve the above equations are as follows. The inlet mass fractions are specified at the channel region (left half of the top edge in Fig. 2). The remaining boundaries are specified as no-flux boundaries for the Maxwell-Stefan equation. For the electron conservation equation, the potential is set to zero at the rib region (right half of the top edge in Fig. 2) and the remaining outer boundaries are defined as no-flux. For the ionic charge conservation equation, the ionic potential is set to ∅ m,0 at the membrane/CL interface and all the other outer boundaries are defined as no-flux.
Physical model for effective specific surface area.-The goal of this study is to geometrically quantify the variation of the effective surface area for varying volume fractions of the CL. Previously, Secanell et al. 15 suggested the following expression for effective surface area by considering the limiting case of zero porosity: where n is the number of agglomerates per unit volume of the CL. However, 15 did not present any physical or geometrical evidence for this expression. Now, n can be expressed as a function of CL porosity and agglomerate dimensions as: Substituting Eq. 23 back in to Eq. 22, it is seen that the effective surface area as modeled by 15 becomes a quadratic function of the catalyst layer porosity because the agglomerate dimensions are assumed to be constant: Equation 24 indicates that a agg is maximized when v = 0.5, and a agg → 0 as v → 0, 1. Apparently, this quadratic dependence of F807 effective surface area on CL porosity was modeled without considering the actual geometrical aspects of how the agglomerates begin to overlap as the CL porosity is decreased, and its consequent effect on agglomerate surface area. Therefore, while Secanell et al.'s expression for a agg 15 is correct for the two limiting cases of v = 0 and 1, there is no physical basis for this expression for 0 < v < 1.
Modeling a agg accurately is of critical importance as it directly affects the agglomerate model predictions for the diffusion loss region. Therefore, a more physically realistic approach based on the geometry of overlapping agglomerates is required to determine how exactly a agg varies with v . This paper describes such an approach based on the work of Taylor et al. 33 who were interested in the effect of biofilm growth on the permeability of porous media. They represented the porous medium as regularly-packed rigid spheres with a uniform film growing on the surface of each spherical particle, and derived expressions for the porosity and specific surface area as a function of film thickness for various sphere-packing arrangements.
By adapting the sphere-packing idea presented in, 33 we developed a new, more accurate method to calculate the effective agglomerate surface area for the classical spherical agglomerate approach. Unlike 33 which determines the porosity as only a function of film thickness, our porosity calculation incorporates the entire input dataset provided for the CL composition.
The catalyst ink composition is specified in terms of the weight fractions of its various components which, consistent with experimental studies, are expressed in terms of the Pt loading. In addition, it should be noted that the CL thickness (t C L ) is assumed to be constant as in. 15 Then the volume fractions of Pt ( Pt ), C ( C ), ionomer ( N ), ionomer in the agglomerate ( agg ), and CL porosity ( v ) are calculated as: where Pt|C and N F P are the Pt and ionomer mass fractions, respectively, in the catalyst powder and the ink. In our model we derive expressions for the reduction in CL porosity and agglomerate surface area due to agglomerate overlap, as opposed to particle-film growth as in. 33 The agglomerate and the ionomer film are assumed to form a spherical particle with radius R = (r agg + δ) which is allowed to overlap with others as shown in Fig. 3. For zero overlap, the effective surface area is obtained trivially as: [30] The porosity value when all the spheres are just touching each other is defined as the critical porosity ( v,cr ) which is specific to each packing scheme. The agglomerates are considered to overlap if v < v,cr . There are four possible arrangements for regularly packed spheres of uniform size: cubic, orthorhombic, tetragonal and rhombohedral, which are named after the shape of the unit cell. 42 For the cubic and orthorhombic packing schemes, the agglomerates start overlapping at porosities of 48% and 40%, respectively, implying that the available surface area, and hence performance, would begin to decline below these porosity values. The corresponding porosities for the tetragonal and rhombohedral schemes are 30% and 26%, respectively. Typical CL porosities are generally less than 50%. Therefore, we preferred to use the tetragonal and rhombohedral schemes (illustrated in Fig. 4) in our study as they provide a higher packing density than the cubic and orthorhombic schemes. Properties for the employed packing schemes   Table III. The packing factor (α m ) in Table III is a dimensionless number that characterizes the packing scheme and is equal to 3 , where V uc is the volume of the unit cell, and D p is the particle diameter.
For v < v,cr , the overlapping ratio ( f ) defines the extent of overlap between particles. It is apparent that increasing f results in decreased void fraction and specific surface area (surface area of spheres/unit cell volume). The variation of porosity and the specific surface area with f depends greatly on the packing scheme as each packing arrangement is characterized by a unique number of contact points between neighboring particles as listed in Table III.
By representing the CL as a bed of regularly-packed spheres, the CL porosity given by Eq. 29 is linked to the void volume fraction of the unit cell. The unit cell's void volume fraction is in turn obtained as a function of the overlapping ratio for the given packing scheme. The unit cell volume for a given packing option is calculated using the packing factor given in Table III, and the void volume is obtained by subtracting the spherical particle volume from the unit cell volume after accounting for the volume lost in the overlapping regions. The total volume loss due to overlapping regions can be expressed in terms of the spherical cap volume and the number of contact points. The where the first term is the volume of spherical particles in the unit cell without accounting for overlapping, and the second term is the total overlapped volume due to the particle's contact with m neighboring particles. Noting that the unit cell volume depends on the overlapping ratio, the CL porosity is obtained from the void volume fraction of the unit cell as: As the CL porosity can be determined from Eq. 29, the goal here is to determine the amount of overlap between the agglomerates for a given CL composition. Equation 32 can be rearranged to yield a third order polynomial for f : Only the smallest root of Eq. 33 satisfies the constraint that 0 < f < 1 while also behaving in a physically realistic manner, and so it represents the desired overlapping ratio. Knowing f , the effective surface area a agg can be derived in a similar manner to Eq. 32. Accounting for the surface area lost due to overlapping particles, the effective surface area is calculated as: [34] It should be noted that Eqs. 32-34 cannot be applied when the CL porosity drops below a limiting value denoted by v,o . For v < v,o , the amount of overlap increases to such an extent that spherical caps from neighboring spheres themselves begin to overlap with each other as illustrated in Fig. 5b. In this case, Eq. 32 underestimates the porosity as it does not account for the overlap of the spherical caps with each other, and Eqs. 32-34 become invalid. A similar situation was examined by 42,43 for Taylor's problem. Adapting their approach, we determined the limiting overlapping porosity v,o for each packing option as listed in Table III. The values of v,o are 6% for the tetragonal-sphenoidal arrangement, and 3% for the rhombohedral arrangement. These porosity values are much too small to be meaningful in the context of a PEM fuel cell CL, and hence it is not worthwhile to solve the relevant higher-order equations for this low-porosity range as in. 42,43 It is readily apparent that a agg = 0 when v = 0, and Eqs. 32-34 allow us to compute a agg for small values of v down to v,o . Hence, for our purposes, it is sufficient to bridge the gap for a agg in the small range 0 < v < v,o by a simple curve fit. A parabolic fit was selected and   [35] In summary, the variation of effective surface area with respect to CL composition is fully described by Eqs. 30, 34, and 35. The variation of effective surface area is plotted against CL porosity in Fig. 6 for the two packing schemes studied here. For each packing scheme, the specific surface area varies linearly with CL porosity for large values of porosity at which the agglomerates do not overlap. The linear behavior is simply due to the number of agglomerates per unit CL volume increasing linearly with the decreasing porosity. The effective surface area reaches its maximum value at the critical porosity value (= v,cr ) when overlap just begins for each packing scheme. Then, due to overlapping agglomerates the effective surface area decreases for v,o < v < v,cr according to Eq. 34. Finally, for ε v < ε v,o , the approximate curve-fit of Eq. 35 is applied. Figure 6 also shows the variation of a agg vs. v as modeled by Eq. 22. 15 Although Eq. 22 is correct at the two limits of the porosity range as stated earlier, the surface area is arbitrarily maximized at 50% porosity and it seriously underpredicts a agg for all other values of ε v compared to the sphere-packing approach presented here.

Validation
The accuracy of our numerical model is validated by comparing its predictions against experimental data. As stated earlier, only the exchange current density and agglomerate dimensions are tuned to match the experimental results of. 31 Utilizing the identical operating conditions and CL architecture of 31 the model (with rhombohedral sphere packing) is shown to accurately reproduce the polarization curve of 31 in the activation, ohmic, and mass transport regimes as shown in Fig. 7.
Additional validation is obtained by comparing the model's predictions (with rhombohedral sphere packing) against experimental data when the CL composition is varied. The effect of Nafion percentage (NFP) on current density at a constant operating voltage as predicted by the model is shown in Fig. 8 along with the experimental data of. 29 As shown in Fig. 8, the agreement is excellent. The current density is seen to increase with NFP up to about 30%, beyond which the decreasing porosity results in a loss of surface area which causes the current density to drop. Figure 8 also shows the calculated porosity value which decreases from 40% to 5% with increasing NFP.
Further validation (with rhombohedral sphere packing) is obtained by examining the current density distribution along the membrane/CL  interface for three different operating regimes. Figure 9 shows that the current density is maximized under the land for the activation region loss (V c = 0.87 V), at the corner of the land for the ohmic loss region (V c = 0.57 V), and under the channel for the diffusion loss region (V c = 0.27 V). This overall trend is in agreement with the results in. 4,13

Results and Discussion
In the previous sections, we have described and validated a new method to calculate the effective agglomerate surface area which plays a critical role in determining the diffusion resistance in the agglom- Figure 9. Current density distribution along the membrane/CL interface for three values of cathode potential, V c (using rhombohedral sphere packing). erate model. Next, we present results from simulations to optimize the CL composition using this new approach to calculate the effective surface area.
The ultimate goal in CL optimization is to achieve high performance at low catalyst loading. For this purpose, one must be able to accurately capture the effect of catalyst loading in numerical models. Hence, we begin by evaluating the performance of our model incorporating the new effective surface area calculation method as the catalyst loading is varied. Figure 10 presents the polarization curves for the different effective surface area calculation methods (the two packing schemes discussed here and Eq. 22) as the Pt loading is varied from 0.2 to 0.4 mg/cm 2 for a fixed Pt|C weight ratio of 20%. Figure 10 indicates that the polarization curve is very sensitive to the effective surface area. The polarization curves in Figs. 10b and 10c obtained with our sphere-packing approach are very similar, except that the rhombohedral scheme predicts higher performances for 0.35 mg/cm 2 Pt loading. In both Figs. 10b and 10c, performance improves substantially as the Pt loading is increased from 0.2 to 0.3 mg/cm 2 . The rate of improvement then slows up to 0.35 mg/cm 2 , beyond which, the performance declines. Therefore, both packing schemes suggest an optimal Pt loading of about 0.35 mg/cm 2 . Figures 10b and 10c match qualitatively the trends seen in experiments. [29][30][31] The drop in performance for Pt loading >0.35 mg/cm 2 in Figs. 10b and 10c can be attributed to decreased porosity at higher Pt loadings which reduces the effective surface area. Such a trend can be observed in the experimental results of 31 for the activation region, and in 29,30 for the diffusion loss region. In contrast, the polarization curves in Fig. 10a obtained with Eq. 22,15 indicate that performance declines monotonically with Pt loading in all the regions of polarization curve. The main reason for this monotonic trend is that Eq. 22 predicts the highest effective surface area at a porosity of 0.5, with the surface area declining continuously for porosities below 0.5 as illustrated in Fig. 6. In our simulation, we have used the CL composition and thickness data provided in; 31 the typical CL porosity is less than 0.5 as reported in. 31 Accordingly, for porosities less than 0.5, a monotonically decreasing performance is observed in Fig. 10a.
To further investigate the behavior of our sphere-packing approach with varying catalyst loading, we examined the effect of Pt|C weight ratio for rhombohedral packing. First, the Pt|C weight ratio is changed from 20 to 40% and the effect of Pt loading on performance is presented in Fig. 11. Figure 11 shows that the performance increases monotonically with Pt loading in contrast to the non-monotonic variation for 20% Pt|C observed in Figs 10b and 10c. Similar to the trend reported in experiments, 29 a higher Pt loading is required to achieve the best performance when the Pt|C ratio is changed from 20 to 40%. Next, the performance of 20 and 40% Pt|C vs. Pt loading are compared at a constant operating voltage of 0.6 V in Fig. 12 as in. 29 In agreement with the experimental trend of, 29 it is seen that Pt|C = 20% consistently performs better than Pt|C = 40% except at very high Pt loading. Therefore, using 20% Pt|C may be a good option for cost-effective catalysts.  Next, we examine the effect of varying the Pt|C weight ratio from 10% to 40% at a constant Pt loading of 0.2 mg/cm 2 in Fig. 13. In agreement with the experimental results of Paganin et al. 30 the performance is seen to vary non-monotonically with Pt|C ratio. In agreement with the results of, 30 smaller limiting currents are observed for high Pt|C ratios. Because the CL thickness is fixed in this simulation, porosity becomes the dominant factor in determining the effective surface area, and hence the limiting current. For Pt|C > 15%, the amount of carbon decreases which increases the overall porosity, which in turn leads to a lower effective surface area, and hence the limiting current is reached earlier.
Next, we investigate the effect of ionomer loading. In this study, the ionomer loading is determined by the ionomer weight fraction (N F P), which is reported in 44 to unify the results in the literature 30,[45][46][47] for ionomer loading. The polarization curves in Fig. 14 are obtained with different effective surface area calculation method as N F P is varied from 25% to 45%. In Fig. 14b, the performance improves with N F P up to 35% for rhombohedral packing scheme, beyond which excessive ionomer content decreases performance by decreasing porosity. Similarly, for the tetragonal packing scheme in Fig. 14c, the optimal value of N F P is 33%. Both packing arrangements suggest optimum N F P values for the CL in agreement with experimental studies. [29][30][31][44][45][46][47][48] Passalaqua et al., 44 Uchida et al. 46 and Russel et al. 31 reported an optimal N F P of 33% which is identical to the value obtained here for the tetragonal packing scheme. In addition, the optimum N F P values deduced from the experimental data of 30 and 45 are between 30-35% and 40%, respectively. In contrast, Fig. 14a shows a monotonically decreasing performance with N F P. This situation can again be explained by porosity values less than 0.5 and the behavior of the effective surface area expression proposed by. 15 The effect of ionomer loading is shown in Fig. 15 for the three effective surface area calculation methods for a fixed operating voltage of 0.6 V. Figure 15 indicates that all three methods predict an optimal value of N F P for the ohmic loss region. The improvement in performance with N F P can be attributed to the improvement in proton conductivity which is proportional to N F P according to the Bruggeman approximation. On the other hand, increased N F P can have an adverse effect on electron transport and the diffusion of reactants. In this study, only the latter effect is modeled, hence the optimal N F P values for the ohmic loss region may be slightly overestimated. Even so, both sphere-packing methods properly capture the performance improvement with increasing ionic conductivity, followed by the adverse effect of decreasing porosity in Fig. 15. The rhombohedral and tetragonal schemes predict optimum N F P values of 35% and 33%, respectively, which are in very good agreement with the experimental results. [29][30][31]44,46 In contrast, the use of Eq. 22 places the optimum N F P at around 27%; moreover, the trend of the curve with Eq. 22 in Fig. 15 also does not match experimental observations. It should be recalled that the main handicap of Eq. 22 is its unphysical   modeling of porosity, whereas the sphere-packing method represents a more realistic approach.
The effect of Pt loading on the optimal value of N F P was also investigated for the rhombohedral packing scheme in Fig. 16 by studying a low (m Pt = 0.25 mg/cm 2 ) and a high (m Pt = 0.5 mg/cm 2 ) catalyst loading case. For the low Pt loading case, the optimal N F P value was 50% as seen in Fig. 16a. For the high Pt loading case in Fig. 16b, the optimal N F P was 20%. In addition, it may be recalled that for the moderate Pt loading of 0.325 mg/cm 2 in Fig. 14b, the optimal N F P was 35%. In general, this suggests that as the Pt loading decreases the optimal N F P value increases, which is consistent with the results reported in. 49 Finally, we investigated the effect of CL thickness by utilizing the rhombohedral packing scheme in Fig. 17. Figure 17 shows that the performance improves with CL thickness from 6 to 10 μm. Beyond 10 μm, the performance falls. As t C L is increased from 6 to 10 μm, the porosity increases from 3.5% to 42% and it promotes reactant transport so the performance improves mainly in the diffusion loss region. On the other hand, charge transfer losses are linearly proportional to t C L . For t C L > 10 μm, the slope of the curves in the ohmic region in Fig. 17 increases, confirming higher charge transfer losses.

Conclusions
The effective agglomerate surface area is a critically important parameter for accurately determining diffusion losses in the classical agglomerate model. A new approach has been presented to model the variation of the effective agglomerate surface area with catalyst layer composition based on two regular sphere-packing schemes: rhombohedral, and tetragonal. It is shown that the sphere-packing methods provide an agglomerate surface area term that varies in a realistic manner with catalyst layer porosity. The results of the new approach are compared with the previously suggested expression in. 15 The effects of Pt and ionomer loading were investigated by using the sphere-packing approach, and the results for both were found to be in excellent agreement with experimentally observed trends. The results confirm that the sphere packing method's success can be attributed to the physically rigorous manner in which the effective surface area is calculated in contrast to the previous expression of. 15 In addition, the model was used to study the effect of Pt|C weight ratio and Pt loading on the optimum ionomer loading, as well as the effect of catalyst layer thickness on performance. In all cases, the model predictions matched experimentally observed trends.

Acknowledgments
This work was conducted under the University of Delaware's Fuel Cell Bus Program to research, build, and demonstrate fuel cell powered hybrid vehicles for transit applications. This program is funded by the Federal Transit Administration at the Center for Fuel Cell Research at the University of Delaware.