Computational Analysis of Gas Transport in Fuel Cell Catalyst Layer under Dry and Partially Saturated Conditions

A computational study is performed to analyze oxygen transport in dry and wet stochastically reconstructed catalyst layers (CLs). CL stochastic reconstructions are generated using random penetrating spheres of a given particle size that agree with the statistical correlation functions of a reference 3D FIB-SEM reconstruction. A nucleation-based full morphology approach is used to partially flood the CLs from within the structure in an attempt to reproduce the conditions in an operating fuel cell. In order to validate the 3D numerical model and reconstruction method, the dry effective diffusivity of CLs and its variation with porosity are obtained and shown to be in agreement with reported literature data. Simulations are then performed for CLs with varying porosity and saturation. Statistical analysis is used to estimate an expression for effective Knudsen radius as a function of porosity and particle size, and then the computed dry effective diffusivities are used to develop a generalized percolation-based correlation function to estimate dry and wet effective diffusivities. The effective diffusivity of CLs with different pore size distribution are also obtained from 3D simulations and using the proposed correlation function and shown to differ by less than 10% for porosities in the range of 0.4-0.7. © The Author(s) 2019. 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.0081907jes]

Water management remains critical to the operation of proton exchange membrane fuel cells (PEMFCs). Low water content in the membrane electrode assembly (MEA) leads to drying of the polymer membrane, thereby leading to a decrease in the protonic conductivity and increase in ohmic losses. 1,2 High water content in the MEA can lead to condensation of liquid water in the channels, gas diffusion layer (GDL), micro-porous layer (MPL) and catalyst layer (CL), causing water accumulation in the porous media. Under these conditions, the gas transport in the MEA is severely hindered and a significant drop in the cell performance is observed. 3,4 Optimizing the water balance to maintain the electrolyte hydrated while removing excess water from the MEA and channels is crucial to improving the performance of the PEMFC at high current densities, 5,6 particularly when the PEMFC is operating at a high relative humidity.
To describe the two-phase characteristics of the PEMFC operation several macro-homogeneous numerical models have been proposed in literature. These models are either based on the solution of a saturation equation derived from the Darcy's law 7,8 or the solution of gas and liquid mass and momentum equations using a closure equation based on the porous media structure such as pore size distribution. [9][10][11] While these models provide valuable insight into the macroscopic effects of the two-phase flow in the PEMFCs, they rely on volume-averaged microstructural information of the porous media and thus, cannot describe some microscale effects, such as loss of pore connectivity on transport and cell performance, as water accumulates in the diffusion media.
Microstructural models can be used to complement the macrohomogeneous models by developing correlations to describe the effective transport properties for the porous media as a function of morphology and local saturation. Several methods, such as pore network model (PNM), [12][13][14][15] lattice Boltzmann method (LBM) [16][17][18][19] and full morphology (FM) or morphological image opening (MIO), [18][19][20][21][22] have been proposed in literature to study the liquid water intrusion in the microstructures of the fuel cell porous media. The most commonly used method is the pore network model (PNM). Pore network models rely on generation of an abstract network of pores and throats with certain radii. These networks can be generated from the real microstructure images obtained using microscopy techniques such as X-ray microcomputed tomography 12 or by using a random network [13][14][15]18 which is calibrated with the experimental data such as porosity profiles, 15 mercury intrusion porosimetry data or saturation-pressure profiles. PNM offers a fast and accurate way to study the effect of localized water effects on the effective transport properties of the porous media. [12][13][14][15] However, generation of the abstract network leads to the loss of the microstructural information such as pore morphology and tortuosity.
Lattice Boltzmann method is a physics based method which has also been employed to study the water transport [16][17][18][19] in porous media. The major advantage of the LBM over the other methods is the ability to identify intricate dynamics of the liquid water and water droplets within the porous media while the other techniques assume a quasistatic water front movement through the porous media. This however, comes at an extremely large computational cost. 19 Full morphology (FM) or morphological image opening (MIO) models describe the water invasion into the porous media, like most PNM, using the Young-Laplace equation. [18][19][20][21][22] The advantage of this technique over the PNM is that it can be used on the real microstructure of the porous media and therefore, further transport and electrochemical simulations can be performed in the partially saturated porous media without any simplifications. It offers a trade-off between the simplicity of PNM and computational cost of LBM. Sabharwal et al. 22 recently showed that this approach can successfully predict water accumulation and effective diffusivities in a GDL sample. Each of the methods described above has been used to study fuel cell porous media such as gas diffusion layer (GDL) [12][13][14][15][16][17][18]23,24 and microporous layer (MPL). 20,[25][26][27] However, there are very limited number of studies on the fuel cell catalyst layers (CLs). Hutzenlaub et al. 28 obtained a CL microstructure using focus ion beam-scanning electron microscopy (FIBSEM) and studied liquid water accumulation on this structure using a PSD based method, where the pores were filled in ascending order corresponding to a hydrophilic (HI) substrate and descending order for the hydrophobic (HO) case. The method of water accumulation proposed in Ref. 28 is questionable, as the advancement of the liquid water front is independent of the liquid capillary pressure. Mukherjee et al. 16 studied the two-phase flow in stochastically generated CL microstructure using LBM and found that at low saturations discrete water clusters are formed which interconnect as the saturation increases. They did not investigate the effect of local saturation on the transport properties and performance. Hannach et al. 29,30 employed PNM to study the two-phase characteristics in the CL. They allowed for water generation from randomly chosen "active agglomerates" in their network and tracked the water interface movement through the porous network and computed the effect of the local saturation on the effective transport properties and performance. Zheng and Kim 31 used LBM to simulate water intrusion based on water production from active ORR sites in a stochastically generated CL microstructure. Zheng and Kim 31 then studied the effective diffusivity of the CL as a function of saturation. However, Zheng et al. 31 only used a single microstructure to study the effective diffusivity at different saturations therefore, statistical variability of microstructures at same and different porosities was not considered. Fathi et al. 32 simulated liquid water intrusion into a stochastic CL reconstruction using a volume of fluid (VOF) approach. Although they computed effective diffusivity for CL reconstructions with different porosities, saturation and domain sizes, the method used in their work to predict the Knudsen effects has previously been shown to underpredict the Knudsen resistance. 31 Thus far, none of the studies have used a FM approach such as the one recently proposed and validated by Sabharwal et al. 22 Further, all of the previous numerical studies 16,[29][30][31][32] used stochastically generated CL structures which were not compared to a real CL microstructure, therefore it is not clear that they are representative of a real CL.
In the present work, an overlapping sphere based algorithm is used to generate stochastic reconstructions of CLs to study the effective diffusivity as a function of porosity and saturation. A method to assess the morphological equivalence of the stochastic reconstructions to a real CL microstructure obtained using focused ion beam-scanning electron microscopy (FIBSEM) is discussed that uses various statistical functions as a measure of its representativeness. Liquid water intrusion is simulated in the CL using a novel nucleation based full morphology algorithm. Voxel based meshes were generated for the dry and wet CL microstructures by extraction of the percolating gas network. The meshes were then used for gas transport simulations that account for both molecular and Knudsen diffusion. Details of the stochastic reconstruction algorithm, water intrusion algorithm, numerical gas diffusion model and statistical functions are provided in Theory section. The results for the comparison of the stochastic reconstructions to FIB-SEM reconstruction and gas transport simulation results are shown in Results and Discussion section. The effective diffusivity was computed for the stochastic reconstructions with different porosities and at different saturations and a correlation for the effective diffusivity is proposed. The simulated values were compared to prior numerical and experimental studies in literature.

Theory
Catalyst layer microstructure reconstruction.-Image based reconstruction.-A thin, low loading CL prepared using inkjet printing, 33 having a platinum loading of 0.025 mg/cm 2 and 30% wt. Nafion loading, was previously imaged using focus ion beam-scanning electron microscopy (FIBSEM) in our earlier work 34 and is used as a reference. The obtained images were used to generate a 3D reconstruction of the CL section with dimensions 848 nm × 447 nm × 1220 nm (424 × 176 × 61 pixels) and a porosity of 39.7%. A detailed description of the image acquisition and image processing algorithms to reconstruct the CL microstructure has been described in our earlier work. 34 Stochastic reconstruction.-Stochastic reconstructions were generated using a random overlapping sphere based algorithm similar to that described in Reference 35. The steps for the overlapping sphere based microstructure reconstructions are as follows: 1. An empty domain ( ) corresponding to a porosity of 100% is generated as a 3D array with user defined dimensions of l × l × l. Building blocks for the solid phase, i.e., spheres of a userdefined radius r d , are generated as local maps using a 3D array A random location is chosen in as a center for map .

3.
map is inserted into the domain by placing its center at the chosen center in the previous step. 4. The process is repeated to choose new centers and insert spheres. 5. As the spheres are continually placed they might start overlapping.
The amount of overlap is controlled in this algorithm using a penetration parameter (ψ) which is used to calculate the region of pixels around the solid sphere in which cannot be used as centers. For 0 ≤ ψ ≤ 1, all pixels within ψr d distance from the surface of the sphere are removed, thereby allowing a penetration of 1-ψ. Therefore, ψ = 0 would allow free penetration of spheres and ψ = 1 would allow no penetration of spheres. 6. Once the desired porosity (ε V ) is reached the algorithm ends. Figure 1 shows a schematic of the different steps described above. As described above, the inputs to the algorithm are the domain size (l), sphere radius (r d ) and porosity (ε V ) of the microstructure. Since the primary objective of this work is to study the gas diffusion in the CL microstructures, the spheres are treated as solid particles without distinguishing between carbon, platinum and ionomer. This is different from the reconstructions in Reference 35 where the spherical particles were treated as carbon particles which were then coated by uniform ionomer films.
Statistical correlation functions.-Two-point correlation function.-The two-point correlation function (S j 2 (r)) is defined as the probability of finding two points at a distance r in the phase j. 36,37 For a digitized porous media, the two-point correlation is computed by traversing a line of length r through the porous media and recording the number of times the ends of the line fall in the same phase. These are then normalized by the total number of traversals. Thus, the two-point correlation function for phase j at r = 0 equals the volume fraction of the phase and the slope of the function is related to the specific interfacial area (interfacial area per unit volume) as follows: 38 , [1] where s j is the specific interphase area, and β is 4 for 2-D, and 6 for 3D images, respectively. 39 Chord length function.-The chord length function (C j (r)) is defined as the probability of finding a line in phase j of length r connecting two interfaces. Therefore, the chord length function provides information about the separation between the phase interfaces. When computed for the void phase in the CL reconstructions, it is representative of the pore size distribution. The chord length function is computed by selecting an interfacial pixel and recording the length of F3067 the chord through the phase j in a given direction to another interfacial pixel. The process is repeated for all interfacial pixels. C j (r) is computed by dividing the number of occurrences of a chord of length r by the total number of occurrences. For this study, the chord lengths are evaluated in the Cartesian directions similar to our earlier work. 34 Chord lengths however, can be evaluated in any direction as shown in Reference 23.
Pore size distribution.-The pore size distribution for the CL microstructures was computed using the sphere fitting algorithm described in Reference 34. To compute the pore sizes, the Euclidean distance transform, i.e., minimum distance of a pore voxel to a solid voxel, is first computed for every pore voxel. The distance transform indicates the radius of the largest sphere that can be placed with the pore voxel as center. The next step is to systematically increase the pore radius r p from the minimum to the maximum value in the distance transform map. For every r p , all pore voxels with distance transform value greater than r p are found. These are treated as the centers for the spheres with neighboring voxels within a distance of r p all being assigned a value of r p . Therefore, every voxel is either the center of a sphere with the radius equal to its distance transform value or part of a larger sphere. Details of the algorithm and implementation can be found in Reference 34.

Algorithm 1
Cluster based full morphology algorithm for water intrusion. 22a Read original input microstructure 0 ; Define the nucleation points, contact angle (θ) and number of steps (n); Compute p c for every pore pixel in 0 using Equation 2 Initialize p in as the minimum liquid pressure to start water intrusion at each nucleation point Initialize l = nucleation points ( l has a value of 1 at the locations where water is present.) Compute p = max(pc )−p in n while s<1 do procedure IDENTIFY LIQUID WATER CLUSTERS Obtain the set of voxels where p c < p in Cluster connected voxels in and store in with each cluster denoted by a unique integer j Identify cluster indices ( j) in connected to l and store in array u end procedure procedure  Record l for mesh generation Increment p in = p in + p end procedure end while a All capital Greek symbols denote 3D arrays with dimensions of input microstructure.
Water intrusion algorithm.-To study the effect of liquid water on the gas transport in the CL microstructure, a quasi-static water intrusion algorithm derived from the full morphology model was developed. 21 The cluster based full morphology (CFM) algorithm used in this study has previously been validated with μ-CT reconstructions of partially saturated GDL reconstructions. 22 However, unlike the GDL where the water was intruded from a single boundary, the liquid water in the CL is intruded using the nucleation mode shown in Figure 2. This is similar to the approach used in the pore network models in References 29,30 where random agglom- erates were activated in the CL to produce water which then intrudes the pores in the CL.
In the CL, the ORR takes place on the Pt surface where water molecules would be produced which might form small clusters. Once these clusters exceed the critical diameter given by the Kelvin equation for the local supersaturation they would continue to grow into bulk liquid. It is expected that such a phenomenon might be occurring in the CL where high local supersaturation might exist especially in the small nano-pores. Therefore, this mode of water injection is deemed to be representative of the water intrusion in the CL. In the present study, pores having 2 nm radius were treated as nucleation sites and provided as input to the algorithm.
The computational implementation for the CFM algorithm is presented in Algorithm 1. 22 A detailed discussion on the algorithm is provided in Reference 22. The key steps are described below: 1. Calculate the pore size distribution for the original microstructure 0 using the sphere fitting algorithm described in Pore size distribution section. 2. Compute the capillary pressure (p c ) required to intrude each pore using the Young-Laplace equation, p c = 2γ cos θ r p [2] where γ is the surface tension of water, θ is the contact angle and r p is the computed local pore radius. For computational implementation, capillary pressure is defined as, where p l and p g are liquid and gas pressures respectively. Therefore, for implementation purposes hydrophilic pores have a negative capillary pressure and hydrophobic pores have a positive capillary pressure. It is assumed that p g is zero. 3. Liquid pressure is then increased using p increments. At each pressure, all pores having capillary pressure less than the liquid pressure are identified as feasible locations where liquid water could exist. 4. A subset of feasible pores is then identified consisting of all pores connected to the liquid water front by clustering the feasible pores and performing a logical AND operation with the liquid water network ( l ). These pores are then intruded with water and the locations are updated in l . 5. The saturation is obtained using, where ε( l ) denotes the volume fraction of voxels with liquid water and ε V denotes the porosity of the microstructure. l is recorded for mesh generation to study the gas transport. 6. Steps 3-5 are repeated until the layer is fully saturated.
The meshes for partially saturated CL reconstructions were obtained by direct conversion of voxels to cells in the mesh. Therefore, the cell size for the mesh was equal to the voxel size. To improve computational efficiency, only the percolating gas network (pores not filled with liquid water) was extracted and used for transport simulations.
Gas transport in CLs.-The gas transport model was described in depth in our earlier study. 34 An overview of the assumptions and key steps is described here.
Oxygen transport in the CL is assumed to be steady state and governed by Fick's law, where D is the overall oxygen diffusivity, c tot is the total gas concentration (assumed constant), and x O 2 is the oxygen molar fraction. Equation 5 assumes that there is no reaction happening in the domain. This assumption is made in order to study the effect of CL geometry on the effective diffusivity. The overall diffusivity, D, is computed by evaluating the Bosanquet equation 40 at every cell, where D molecular is the molecular oxygen diffusivity and D Kn is the Knudsen diffusivity. The Knudsen diffusion coefficient is calculated using, 35,41,42 where r p is the local pore radius obtained using the sphere fitted pore size distribution 34 and therefore, changes in every cell, R is the universal gas constant, T is the absolute temperature of the gas and M O 2 is the molecular weight of oxygen. For pores filled with liquid water the diffusion coefficient was assumed to be zero because of the four orders of magnitude difference between the molecular oxygen diffusion coefficient in air and water. 43 The boundary conditions used to simulate the oxygen transport in the microstructure are, and (Dc tot ∇x O 2 ) · n = 0 everywhere else, [8] where 1 is the inlet plane and 2 is the outlet plane opposite to the inlet plane. The effective diffusion coefficient (D e f f ) of the porous media was computed using whereṄ out O 2 is the total molar flow rate of oxygen at the outlet plane 2 , L is the shortest distance separating inlet and outlet planes, A is the cross-section area and, c in O 2 and c out O 2 are the concentrations of oxygen at the inlet and outlet boundary faces, respectively.
The effective diffusion coefficient calculated from Equation 9 ignores the reaction in porous media. Several volume averaging studies have shown that the diffusion coefficient can be affected by the heterogeneous reaction. [44][45][46] In order to account for the reaction, a volumeaverage model such as the one proposed by Whitaker 45 for first-order heterogeneous reactions in porous media would have to be developed. For heterogeneous reactions, as is the case for fuel cell CLs, Whitaker 45 showed that, assuming a first-order, constant reaction rate reaction, the reaction term leads to a similar equation, with an effective diffusion coefficient that is purely geometrical in nature, similar to the value estimated in this article, and an additional term, due to the reaction, which manifests itself as a convective transport term in the volume averaging. 44,45 The latter term is however challenging to evaluate in fuel cells. Whitaker's analysis 45 assumed a constant effective reaction rate, however this is not the case for the electrochemical reaction in fuel cells as the reaction rate is a function of the overpotential. Additional work is therefore required to derive a volume averaged expression to express reactive gas transport in fuel cell CLs as the expressions used in the porous media community 44-49 cannot be used directly. For first and second-order reactions, the convective term has been shown to be significant only at high Thiele modulus. [44][45][46]

Results and Discussion
Stochastic reconstructions.-The overlapping sphere based algorithm described in Stochastic reconstruction section was used to generate multiple stochastic reconstructions to study the gas diffusivity as a function of porosity and saturation. However, it is important to determine the statistical equivalence of the stochastic reconstructions compared to the structure of a real CL. Statistical correlation functions such as pore size distribution, two-point correlation function in the void phase (S v 2 (r)) and chord length function in the void phase (C v (r)) were used as parameters to quantify the representativeness of the stochastic reconstructions compared to the real CL microstructure obtained from FIBSEM.
For the comparison, stochastic reconstructions with a porosity of 40% similar to the porosity of 39.7% from the FIBSEM microstructure were used. A domain size of 600 nm × 600 nm × 600 nm with a voxel resolution of 2 nm in every direction was used. The domain size was calculated based on the findings from our previous work 34 where a domain size of nearly 500 nm in every direction was found to provide a representative elementary volume for gas transport. Since the particle radius for the spherical particles used to generate the stochastic reconstructions was not known, a parametric study was performed on the spherical particle radius to determine the value that would better represent the real CL structure. The spherical particle radius was varied between 20-50 nm (10-25 voxels), based on the reported size of a primary carbon particle 50 and 10 reconstructions were generated for each set of parameters to ensure statistical significance of the results. In the current study, the particle size for each reconstruction was considered constant even though in reality the carbon particles would have a particle size distribution. 50 The effect of multiple particle sizes will be evaluated in a future study. For this study, the penetration parameter (ψ) was set to zero allowing free penetration of the spheres so that only the particle radius was the unknown parameter. Although not shown here, limiting the degree of penetration, i.e., reducing the amount of overlap permitted between spheres, results in a decrease in the pore sizes. Statistical functions were computed for each of the reconstructions with different r d values and compared to those obtained from the FIBSEM dataset.
Figures 3a-3d shows a stochastic reconstruction with 40% porosity using particles with a radius of 20 nm, 30 nm, 40 nm and 50 nm respectively as well as the local pore size. The local pore radius is computed using the pore size distribution algorithm described in Pore size distribution section and is used to account for Knudsen effects in the gas diffusion simulations. Figure 4 shows the pore size distribution (PSD) for the stochastic reconstructions. Since Figure 4 shows the pore size distribution as a probability distribution function, the area under the curve for the four figures is 1. For all the cases, the smallest pore size is determined by the voxel resolution, i.e., 2 nm. The largest pore radius increases with an increase in the particle size. For the reconstructions shown in Figure 3, the largest pore radius increases from 38 nm to 55, 67 and 90 nm as the particle radius is increased from 20 nm to 50 nm. The increase in pore size becomes more evident when analyzing the pore size distribution for the 10 different reconstructions for each of the particle sizes. Further, the probability of finding pores larger than 100 nm in diameter is 0 for particle radius of 20 nm but increases to 1%, 7% and 20% for increasing particle radius.
The red dashed line in Figure 4 shows the pore size distribution of the reference CL reconstruction obtained using FIBSEM. It can be seen from Figure 4 that the stochastic reconstructions with particle sizes of 20 nm and 30 nm overpredict the fraction of small pore sizes compared to the FIBSEM reconstruction. Stochastic reconstructions with particle size 40 nm show similar pore size distributions to that of   . Two-point correlation in the void phase in: a) x and b) y direction for stochastic reconstructions with particle radius of 20 nm; c) x and d) y direction for stochastic reconstructions with particle radius of 30 nm; e) x and f) y direction for stochastic reconstructions with particle radius of 40 nm; g) x and h) y direction for stochastic reconstructions with particle radius of 50 nm. the FIBSEM reconstructions. They even show some of the extremely large pores with diameter greater than 150 nm. With a further increase in the particle radius to 50 nm, the pore sizes become much larger than those observed in the FIBSEM reconstructions with some pores becoming as large as 230 nm in diameter. Figure 5 shows the two-point correlation function in the void phase (S v 2 (r)) in the x and y directions for stochastic reconstructions with spherical particle radius of 20 nm to 50 nm. The red dashed line in Figures 5a-5h shows the two-point correlation for the FIBSEM mi-crostructure. The two-point correlations in the x and y directions for reconstructions with particle size of 40 nm are similar to those obtained from the FIBSEM microstructure in the x and y direction. For reconstructions with particle radius of 20 nm and 30 nm, the two-point correlation function is underpredicted for r < 100 nm. For reconstructions with particle radius of 50 nm, the two-point correlation shows slightly higher values in the x direction but similar values in the y direction. It must be noted that the two-point correlation function in the z direction is not used for comparison because of the artificial anisotropy of the FIBSEM microstructure in the FIB slicing direction. Due to a lower resolution of nearly 20 nm in the FIB slicing direction (z) there is loss of information in the z direction of the CL microstructure which results in the artificial anisotropy. This has been discussed in detail in our previous work. 34 As discussed in Two-point correlation function section, the slope of the two-point correlation function at r = 0 is proportional to the interfacial area between the solid and pores. From Figures 5a-5d, it can be seen that the slopes of the stochastic reconstructions with particle radius of 20 nm and 30 nm are steeper than the FIBSEM reconstruction. This indicates that the specific interfacial area between the solid and void phase for these reconstructions would be higher than that observed in the FIBSEM reconstruction. However, the slope of the two-point correlation function for reconstructions with particle radius of 40 nm are very similar to those from the FIBSEM reconstructions as seen from Figures 5e and 5f. The slope of the two-point correlation function for reconstructions with particle radius of 50 nm appears lower than the slope of the FIBSEM. The average specific interfacial area for the 10 reconstructions for each particle radius are shown in Table I. The stochastic reconstructions with particle radius of 20 nm and 30 nm overestimate the total specific interfacial area by 88% and 24% respectively. This can be explained by the higher fraction of smaller pores in these reconstructions compared to the FIBSEM reconstruction. The average total specific interfacial area for stochastic reconstructions with particle radius of 40 nm is much closer to the FIBSEM reconstruction with a difference of only 7%. With a further increase in the particle radius to 50 nm, the average total specific interfacial area decreases further due to increase in the pore sizes as shown in Figure 4. Figure 6 shows the chord length function in the void phase (C v (r)) for stochastic reconstructions with particle radius of 20 to 50 nm. The chord length function for stochastic reconstructions with particle radius of 40 and 50 nm closely resembles the functions from the FIB-SEM data. For stochastic reconstructions with particle radius of 20 and 30 nm, the chord length function has a higher value for smaller r compared to the FIBSEM data. Table I shows the average mean chord length (λ) for the reconstructions calculated using, [10] The mean chord length for the reconstructions with particle radius of 40 nm are nearly identical to those from the FIBSEM whereas the mean chord lengths for the reconstructions with r d of 20 and 30 nm are significantly smaller. Although the chord length function appears to be similar between the reconstructions with particle radius of 50 nm and the FIBSEM reconstruction, the mean chord length for this case is 65.6 nm which is nearly 20% higher than the mean chord length of the FIBSEM data. This is likely due to a higher probability of the longer chords due to the larger pore sizes as seen in Figure 4d. The correlation functions and pore size distributions for the different reconstructions with the same particle radius are also similar thereby, indicating that the reconstructions can be considered as REVs.
Comparison of the statistical functions of the stochastic reconstructions to those obtained from the FIBSEM reconstructions shows that: a) selecting the appropriate particle size in reconstructions is critical to achieve a statistically equivalent structure to a real CL; and b) in this case, the stochastic reconstructions with particle radius of 40 nm are statistically equivalent to the FIBSEM microstructure of the CL used in this study. Therefore, reconstructions with a particle size of 40 nm are used for analysis of diffusion as a function of porosity and local saturation in the CL.
CLs gas diffusivity under dry conditions.-Effect of porosity.-To determine a correlation for the effective diffusivity of CLs, the effect of porosity on the effective diffusivity is analyzed. The overlapping sphere based algorithm described in Stochastic reconstruction section was used to generate multiple reconstructions with varying porosities and r d = 40 nm. The domain size for each reconstruction was 600 nm × 600 nm × 600 nm with a voxel resolution of 2 nm in every direction. The porosity of the reconstructions was varied between 1-10% in increments of 1% and from 10-100% in increments of 10%. 10 reconstructions were generated for each porosity value to provide statistically significant results. To simulate gas diffusion in the generated stochastic reconstructions the percolating pore network is extracted.
Reconstructions with porosities below 6% did not contain a percolating pore network in at least one or more Cartesian directions. Therefore, a porosity of 5% is chosen as the percolation threshold (ε th ), which is the porosity above which a percolating pore network exists. Similarly, it must be noted that the stochastic reconstructions with porosity higher than 70% were found to have no percolating solid network thereby, indicating that such a structure might be infeasible. Therefore, the effective diffusivity was only computed for stochastic reconstructions with porosities upto 70%. Figure 7 shows the pore phase for one of the CL reconstructions with porosity of 7%, 50% and 70%. An increase in porosity leads to an increase in the pore sizes, with the largest pore radius increasing from 35 nm at a porosity of 7% to 92 nm at a porosity of 50% and 108 nm at a porosity of 70%. Figure 8 shows the effective diffusivity in the x, y and z directions as a function of porosity for the different stochastic reconstructions. The effective diffusivities are computed at a temperature of 80 • C using molecular oxygen diffusivity of 0.273 cm 2 /s. As expected, the effective diffusivity in the Cartesian directions increases with an increase in porosity.
It can be seen from Figure 8 that for a given porosity the values of the effective diffusivities in the three directions are nearly identical indicating an isotropic structure. This is expected because the overlapping sphere algorithm is not biased in any direction and the building blocks are spheres which are also isotropic. Further, it can also be seen that the effective diffusivities for the 10 reconstructions with the same porosity show similar effective diffusivities with a standard deviation of less than 3.5% about the mean value. This also ensures that the average diffusivities computed for each porosity using 10 reconstructions is statistically representative.
The results from Figure 8 can be used to develop a functional dependence of the effective diffusivity of the CL to its porosity. Macro-scale MEA models 9-11,51,52 rely on such semi-empirical relations to compute the effective diffusivity in the fuel cell porous media. Recently, Zheng and Kim 31 and Shin et al. 53 have proposed a tortuosity model given by Equation 11 to estimate the effective diffusivities for the CLs, where the tortuosity factor (τ) in Equation 11 is related to the porosity using a Bruggeman type correlation, where α is a fitting parameter which has a value of 0.5 for the Bruggeman correlation. Zheng and Kim 31 computed a value of 0.75 for α and used D = 1.07D molecular to fit the effective diffusivity values. Shin et al. 53 computed the value of α to be 0.746 using the shape of streaklines and used D = D Kn avg . However, a common problem with the correlation of the form proposed in Equation 11 is the lack of information about the connectivity of the pores, as the porous media will F3072 Journal of The Electrochemical Society, 166 (7) F3065-F3080 (2019) Figure 6. Chord length function in the void phase in: a) x and b) y direction for stochastic reconstructions with particle radius of 20 nm; c) x and d) y direction for stochastic reconstructions with particle radius of 30 nm; e) x and f) y direction for stochastic reconstructions with particle radius of 40 nm; g) x and h) y direction for stochastic reconstructions with particle radius of 50 nm. continue to have a non-zero diffusivity even at extremely low porosities. A consequence of this can be found in the optimization study by Secanell et al. 51 where the optimal CL porosity was computed as 2.5%, whereas in reality such low porosities are likely to lead to a loss of percolating pore network.
An alternate correlation to determine the effective diffusivity for the CLs is based on percolation theory 11,54-56 and given by, ) unless CC License in place (see abstract  where ε th is the percolation threshold value and μ is a fitting parameter. The Heaviside step function (H(ε V −ε th )) in Equation 13 is added to enforce that, for porosities below the percolation threshold, the effective diffusivity is zero because a percolating gas diffusion network would not exist. This is a major distinguishing factor from the Bruggeman correlation and previous correlation functions based on Equation 11 such as the work of Zheng and Kim 31 and Shin et al. 53 which could predict non-zero effective diffusivities below the percolation threshold.
The bulk diffusivity (D pore ) of the CL in Equation 13 is obtained using the Bosanquet approximation given by Equation 6 where the average Knudsen diffusion coefficient (D Kn avg ) can be computed using Equation 7 and the local pore radius (r p ) must be replaced by a mean pore radius (r avg ). Equations 13 and 6 are fitted to the effective diffusivities of the stochastic CL reconstructions shown in Figure 8 to estimate exponent (μ).
The statistical reconstruction simulation results can also be used to estimate the fitting parameters in Equation 13 and the most appropriate Knudsen mean pore radius. In this case, ε th is estimated by direct analysis of the percolating pore network for the stochastic reconstructions. As mentioned earlier, for porosities below 6% there was loss of percolating pore network in at least one or more Cartesian directions for the reconstructions. Therefore, ε th is chosen to be 0.05 which is the value above which a percolating pore network exists.
To compute the average Knudsen diffusion coefficient (D Kn avg ), an average pore radius (r avg ) must be computed. The average pore radius is computed as the arithmetic mean of the pore radii for every pore, 31,35,53 r avg = i r i dX i dr dr [14] where dX i is the volume fraction of pores having radius r i . Since the pore radius is directly proportional to the Knudsen diffusivity,  this method assumes that the average Knudsen radius is such that it results in an average Knudsen diffusivity equivalent to the effective Knudsen diffusivity obtained when transport through all pores occurs in parallel. Figure 9 shows the average Knudsen radius computed using the PSD of the stochastic reconstructions with r d = 20, 30 and 40 nm and Equation 14. As discussed earlier, for porosities below 6% and higher than 70% there is a loss of percolating pore and solid network respectively. Therefore, r avg for these porosities are not shown in Figure 9 and not taken into account to determine the best fit curve. The best fit curves in Figure 9 were estimated by minimizing the leastsquare of the difference between r avg . r avg = r d 1.66ε 1.65 V + 0.289 [15] For the coefficients given in Equation 15, the R 2 value for the fit is 0.996. Once ε th and r avg have been obtained, the value of μ is obtained by minimizing the least-square difference between the results in Figure 8 and Equation 13 with ε th = 0.05. The optimal value of μ obtained for all the cases is 1.90. The R 2 values for the fit is 0.994 indicating a good fit between the data and Equation 13 with ε th = 0.05 and μ = 1.90. In summary, Equation 13 with ε th = 0.05, μ = 1.90 and D pore obtained using Equations 6, 7 and 15 provides an excellent prediction of gas diffusivity for a statistically representative CL structure. Figures 8a-8c show the best fit for effective diffusivity in the x, y and z directions as a function of porosity. Figure 8d shows the best fit line for the average effective diffusivity (D e f f ) of the reconstructions with r d = 40 nm as a function of porosity.
10 reconstructions with r d = 20 nm and porosities ranging from 0.06-0.7 were generated. The effective diffusivity for these reconstructions was computed as shown in Figure 8d. The effective diffusivity for reconstructions with r d = 20 nm are smaller than the reconstructions with r d = 40 nm due to the smaller pore sizes for the former case. This was also shown earlier in Figure 4 for reconstructions with a porosity of 0.4. To further assess the validity of the Equations 6, 7 and 15, they were used to predict the effective diffusivity for reconstructions with r d = 20. From Figure 8d it is clear that the correlation proposed by Equations 6, 7 and 15 with μ = 1.90 and ε th = 0.05 can accurately predict the dry effective diffusivity even for reconstructions with a different microstructure from those used to develop the correlations.
Comparison to literature.-In order to assess the validity of the diffusivities calculated using the stochastic reconstructions and used to derive the correlation function above, the predicted diffusivities in this article are compared to previously reported literature values. To aid in this comparison the formation factor (F ) is used. The formation factor is defined as, [16]  35 Fathi et al., 32 Siddique and Liu, 57 Zheng and Kim 31 and Shin et al. 53 computed the effective diffusivities numerically using stochastic reconstructions for CLs. Yu et al. 58 and Inoue et al. 59 experimentally measured the effective diffusivity of CLs with different ionomer loadings.
The average formation factor value at a given porosity is determined by averaging the effective diffusivities in the three Cartesian directions over the 10 reconstructions. The average formation factor for the stochastic reconstructions at a porosity of 40% (similar to the porosity of 39.7% for the FIBSEM CL) is 0.037 which is similar to the average formation factor of 0.025 ± 0.010 (x and y average and standard deviation) for the reference FIBSEM CL. The standard deviation in the formation factor for the reference FIBSEM CL is computed from the domain size study performed in Reference 34. The z direction values for the FIBSEM reconstruction are not considered due to the artificial anisotropy in the z direction. 34 Figure 10 provides a comparison of the obtained average overall formation factors obtained from the stochastic reconstructions to those previously reported in literature. 31,34,35,53,[57][58][59] Yu et al. 58 measured the effective diffusivities for CLs with platinum supported on amorphous carbon and platinum supported on graphitized carbon with different ionomer to carbon ratios (I/C) ranging from 0.5-1.5. The porosity and pore size distribution was measured using mercury intrusion porosimetry (MIP). The formation factors obtained from the stochastic reconstructions show good agreement with the predicted effective diffusivities with an error of less than 25%. Figures 11a-11c show the pore size distributions of the stochastic reconstructions with porosities of 0.3-0.5 compared to those measured by Yu et al. 58 for similar porosities. For all porosities, the PSDs from Yu et al. 58 are shifted toward higher pore sizes compared to the stochastic reconstructions. Further, the PSDs from Yu et al. 58 did not consider pore sizes less than 25 nm due to limitations of the equipment. This resulted in the higher probability values for the larger pores compared to the PSDs from the stochastic reconstructions. The discrepancy in the formation factor values is likely due to the different microstructures for the CLs from the two studies as evident from the PSDs.
Inoue et al. 59 measured the effective diffusivity of CL samples with ionomer to carbon ratios of 0.8-1.4. The formation factor values reported by Inoue et al. 59 are 3-4 times smaller than predicted effective diffusivity values from the stochastic reconstructions with porosities in the range of 0.4-0.6. For porosities higher than 0.6, the error reduces to less than 40%. To identify the reason for discrepancy in the formation factor values, the pore size distributions reported by Inoue et al. 59 using N 2 adsorption are analyzed for ionomer to carbon ratio of 1 and 1.3 corresponding to a porosity of 0.5 and 0.6 respectively. Figures 11c  and 11d show the comparison of the PSDs from Inoue et al. 59 to the stochastic reconstructions. For ε V = 0.5, the PSD from Inoue et al. 59 are similar to the stochastic reconstructions with pore sizes as high as 213 nm. For ε V = 0.6, the PSD from Inoue et al. 59 shows a higher probability for pore sizes less than 100 nm compared to the stochastic reconstructions, which might lead to a lower effective diffusivity observed by Inoue et al. 59 Analysis of the two PSDs however, does not provide any conclusive information to explain the high discrepancy between the predicted formation factors in this study and those obtained by Inoue et al. 59 Further, the formation factor calculations on FIBSEM reconstructions of the CLs performed by Inoue et al. 59 also overpredicted the formation factors by 2 times. Since Inoue et al. 59 did not report any of the statistical functions for their reconstructions it is difficult to ascertain the reason for the overprediction of the diffusivity values by the simulations.
Lange et al. 35 numerically studied the effective diffusivity in stochastic reconstructions generated using overlapping spheres similar to this work. The formation factor values from Lange et al. 35 are much higher than those obtained in this work. This is surprising because the particle diameter used by Lange et al. 35 is 20 nm which is smaller than the one used in this study and the effective diffusivity increases with an increase in particle diameter. 35 Further, Lange et al. 35 also accounted for local Knudsen effects similar to this work. The difference between the effective diffusivity values could be attributed to the method used to compute the local Knudsen effects. Lange et al. 35 used an average length in the directions around a pore. This approach underestimates the Knudsen resistance as shown by the Zheng and Kim 31 who compared the approach proposed by Lange et al. 35 to results from erosion-dilation based pore radius (similar to this study) and LBM.
Fathi et al. 32 used an overlapping sphere based stochastic reconstruction, similar to Lange et al. 35 and this article, to study the effective diffusivity as a function of porosity. As seen in Figure 10, the diffusivity values obtained by Fathi et al. 32 are higher than those obtained in the present study. This is again due to the underprediction of Knudsen effects due to the use of an average length around the pore similar to Lange et al. 35 Siddique and Liu 57 used a heuristic pixel based reconstruction method to generate CL reconstructions with different porosities and then performed numerical simulations on the reconstructions to obtain the effective diffusivities. The results from this study are in good agreement with those from Siddique and Liu 57 as shown in Figure 10 with a discrepancy of less than 10% for porosities in the range of 0.4-0.6.
Zheng and Kim 31 numerically computed the effective diffusivity for CL microstructures with different porosities generated using a sphere based simulated annealing method. The effective diffusivities by Zheng and Kim 31 were computed using continuum simulations using a local Knudsen resistance with the pore radius computed using morphological image opening, similar to the spherical pore radius used in this study. However, for the continuum simulations they used a mean pore radius to estimate the Knudsen diffusivity and did not account for the local Knudsen resistances. They also used LBM to compute the effective diffusivities of the microstructures. The average formation factor values from the stochastic reconstructions at different porosities in this study are in good agreement with the formation factor results from Zheng and Kim 31 using LBM and erosion-dilation approach with a maximum error less than 10% except for the formation factor computed at a porosity of 30% using the erosion-dilation method where the error is 22%. This comparison also shows the accuracy of continuum models that account for local Knudsen effects is similar to that of LBM while being computationally efficient.
Shin et al. 53 computed the effective diffusion coefficient for randomly generated CL microstructures. Although the formation factor values from Shin et al. 53 are in good agreement with the results from Zheng and Kim 31 and this work, the numerical method used by Shin et al. 53 is questionable. Firstly, the voxel resolution was set to 20 nm which would ignore the smaller pores in the CL which are seen in the FIBSEM reconstruction 34 as well as experimental PSDs. 58,60,61 Coarsening the voxel resolution also leads to an increase in the effective diffusivity value as shown in our previous work. 34 Secondly, the  [17] where V is the porosity, τ is the tortuosity factor computed from an effective relation between the effective molecular diffusivity to the bulk molecular diffusivity 53 and D Kn avg is the average Knudsen diffusivity computed from a mean pore size. This correlation ignores any effect of the molecular diffusivity which is typically accounted for using the Bosanquet equation. 31,34,35 Further, local Knudsen effects are ignored when using a mean pore diameter which would lead to higher diffusivities. 31,35 This is partially compensated by multiplying the average Knudsen diffusivity by the ratio of ε V τ . Comparison of the average formation factors from the stochastic reconstructions used in this work to previously reported numerical and experimental data shows that the values are in good agreement. Additionally, the previous numerical studies did not compare the statistical equivalence of the generated microstructures to an actual CL morphology which was done for the stochastic reconstructions used in this study. This provides good confidence in the predicted diffusivity values to develop an effective diffusivity correlation for the CL.
CLs gas diffusivity under wet conditions.-To study the effect of local saturation on the effective diffusivity, liquid water was intruded into the reconstructions using the nucleation based water intrusion method described in Water intrusion algorithm section. The smallest pores, i.e., pores having a radius of 2 nm, were treated as nucleation sites and assumed to be initially filled. The reason for treating these pores as nucleation sites has been explained earlier in Water intrusion algorithm section. Water intrusion was then carried out using the CFM algorithm described in Algorithm 1 with n = 100 and a hydrophobic contact angle (θ) of 93 • was assumed based on the environment scanning electron microscopy measurements reported in Reference 62. Although a contact angle of 93 • was assumed for this study, it was shown in our previous work 22 that the liquid water distribution using the CFM algorithm does not change with contact angle as long as the pore wettability (i.e., hydrophilicity or hydrophobicity) remains the same. The contact angle only changes the value of the capillary pressure required to intrude the pores. The liquid distributions, and hence the effective diffusivity at a given saturation, does not change with contact angle. Partially saturated CL reconstructions were recorded at intervals of at least 5% saturation difference. The partially saturated CL reconstructions were then converted to voxel based meshes by direct conversion of voxels to mesh cells. The pores filled with liquid water were assumed to be impervious to gas transport and therefore, considered part of the solid region. The local pore radii were recomputed to account for the change in the Knudsen effects due to liquid water intrusion. The effective pore networks, i.e., pores not flooded with liquid water, were extracted and used for gas transport simulations.
Effect of local saturation.-To study the effect of local saturation on the effective diffusivity, stochastic reconstructions with r d = 40 nm and porosities in the range of 0.3-0.6, typical for fuel cell CLs, 58,63 were used. Figure 12 shows the liquid water distribution in one of the CL reconstructions with a dry porosity of 0.5 at various saturations. Figure 13 shows the ratio of average wet effective diffusivity to the average dry effective diffusivity for stochastic reconstructions with porosities in the range of 0.3-0.6. As shown earlier for the dry case, the effective diffusivities in the three Cartesian directions were similar therefore, only the average effective diffusivity is shown for the wet case. The wet diffusivity of the CLs decreases with an increase in the saturation as seen in Figure 13. Further it can be seen that for a given saturation the ratio of the wet effective diffusivity to the dry effective diffusivity is almost identical for all the stochastic reconstructions irrespective of their porosity. This indicates that the ratio of the wet effective diffusivity to the dry effective diffusivity is independent of the dry porosity. The predicted effective diffusivities of the partially saturated stochastic CL reconstructions are in good agreement with the average (of x and y) wet diffusivities obtained from the FIBSEM CL reconstruction with a maximum error of 18% for saturations below 45%. Similar to the dry case, the results in Figure 13 can be used to estimate a correlation for the effective diffusivity of partially saturated CLs. A power law is commonly used to describe the ratio of wet effective diffusivity to dry effective diffusivity of fuel cell porous media. 23,31,32,64 Equation 18 shows the correlation used to estimate the best fit curve in Figure 13, where D e f f dry is given by Equation 13, s is the saturation and γ is an empirical coefficient estimated by minimizing the least square difference between the results in Figure 13 and Equation 18, The optimal value of γ estimated using a non-linear solver is 2.41 which gives a R 2 value of 0.995 for the curve fit. To estimate the functional form of the Heaviside step function in Equation 18, the saturation threshold (s th ), i.e., the saturation above which a percolating gas network does not exist, is shown in Figure 14a. The saturation threshold increases from 0.84 to 0.92 when the porosity increases from 0.3-0.6. Thus, a Heaviside step function of the form H(s − s th ) cannot be used as s th depends on the porosity of the CL. However, a plot of the effective porosity at the saturation threshold shown in Figure 14b, is nearly constant at 0.047 except for the reconstructions with ε V = 0.5, which is slightly higher at 0.056 due to outliers. Therefore, the functional form of the Heaviside step function shown in Equation 18 is appropriate as it accounts for the loss of the percolating gas network with an increase in saturation with the percolation threshold (ε th ) still having a value of nearly 0.05. Figure 13 also shows a comparison of the predicted wet effective diffusivities from the present study to previous numerical studies. 31,32 Fathi et al. 32 intruded water into the stochastic CL reconstructions using a nucleation approach similar to this study. However, they used a physics based volume of fluid (VOF) method to track the liquid intrusion as opposed to the image based CFM approach used in this study. Fathi et al. 32 proposed a correlation for the ratio of wet effective diffusivity to dry effective diffusivity of the form shown in Figure 13. The wet effective diffusivities obtained by Fathi et al. 32 are higher than those predicted from the current study. The higher values obtained by Fathi et al. 32 are likely due to the underprediction of the Knudsen effect as shown earlier for their dry effective diffusivities. Further, Fathi et al. 32 did not account for the change in Knudsen radius with saturation which results in an overprediction of the effective diffusivities as shown by Zheng and Kim. 31 Zheng and Kim 31 used LBM to intrude liquid water into a stochastic CL reconstruction. They claimed to have simulated liquid water generation due to the ORR in the CL however, this is difficult to assess as the numerical details for the liquid water intrusion model were not provided. The wet effective diffusivities computed using LBM by Zheng and Kim 31 and fitted to a power law, similar to this study, are shown in Figure 13. The value of γ computed by Zheng and Kim 31 was 3.062 for a pressure of 1 atm and changed to 2.982 at a pressure of 1.5 atm. Zheng and Kim 31 also reported that when using pore scale simulations without accounting for local Knudsen effects the γ value reduced to 1.815. The value of γ obtained in the current study is smaller than that predicted from the LBM simulations in Reference 31.
The discrepancy between the γ value from Zheng and Kim 31 and the current study might be due to two reasons. Firstly, the higher order LBM method is able to account for Knudsen effects more accurately for water-air interfaces than the assumption of impermeable water interface used in this study. The second reason could be the difference in the microstructure of stochastic reconstructions. The mean pore diameter for the reconstructions by Zheng and Kim 31 ranges from 71-154 nm for porosity range of 0.3-0.7 while for the present study the mean pore diameter ranges from 50-100 nm for the same porosity range. Statistical correlation functions were not provided by Zheng and Kim 31 for their structures therefore, the differences in the morphology of the reconstructions cannot be analyzed. Further, Zheng and Kim 31 used only a single base structure to compute the wet effective diffusivities so they did not account for any statistical variability in the structures or effect of porosity. In the present study, the correlation is derived based on nearly 450 data points from 40 reconstructions with porosities in the range of 0.3-0.6. The value of γ computed from the present results can range between 2.24-2.51 depending on the data set used. Therefore, it is reasonable to assume that the wet effective diffusivities predicted in the current study are within the bounds of microstructure variation.
Thus far, the article has shown that using an effective porosity for the wet samples defined as ε e f f = ε V (1 − s), resulted in a constant saturation threshold. It is reasonable to ask if, instead of using Equation 18, which adds a new term and coefficient to Equation 13, the latter could be used by simply replacing ε V by ε e f f in Equations 13 and 15. Thus, the following correlation, is used with μ = 1.90 and ε th = 0.05 to predict wet diffusivity. Figure 15 shows a comparison of the two correlations given by Equations 18 and 19 as a function of the effective porosity. For Equation 18, the wet effective diffusivity depends on the dry effective diffusivity and hence, the dry porosity therefore, the two extreme cases with porosities of 0.3 and 0.6 are shown in Figure 15. The expression given by Equation 19 provides a R 2 value of 0.986 compared to 0.995 when using Equation 18. Therefore, Equation 19, although does not provide as good of a fit as Equation 18, provides an excellent approximation while removing the need for an extra parameter. The correlation for the effective diffusivity proposed in this study was developed by validating the statistical functions of the stochastic reconstructions against an inkjet printed CL. Further, the predictions of the dry diffusivity of the CLs with different porosities were found to be consistent with experimental data which were for different CLs than the one used in this study. Also, based on comparison of the predictions for the reconstructions with a different r d (than the one used to develop the correlation), as shown in Figure 8d, the correlation should be valid for CLs with different particle sizes (r d ). In summary, Equation 19 combined with Equations 6, 7 and 15, with only two fitting parameters, provides an accurate prediction of dry and wet diffusivities in CL and is suitable for use in volume-average single cell models.

Liquid-vapor interfacial area in the CL.-
The liquid-vapor interfacial area variation with saturation is an important parameter used to model the liquid water evaporation in the macro-homogeneous models. 11,65,66 Zhou et al. 65,66 showed that under wet conditions evaporation of the liquid water in the CL becomes crucial to transport liquid water out of the CL and MPL and prevent complete flooding. Therefore, knowledge of the liquid-vapor interfacial area variation with saturation from the micro-scale models can be used in the macrohomogeneous models to accurately describe the evaporation process.
As described earlier in Effect of local saturation section, stochastic reconstructions with r d = 40 nm and porosities ranging from 0.3 to 0.6 at different saturations were used. Figure 16 shows the specific liquid-vapor interfacial area as a function of the saturation. The specific interfacial area is defined as the liquid-vapor interfacial area per unit volume of the CL. Since the total CL volume for all the reconstructions irrespective of their porosity is the same, the specific interfacial area is directly proportional to liquid-vapor interfacial area. The interfacial area increases with saturation to a maximum at a saturation of nearly 75% and decreases to zero at 100% saturation. An increase in porosity leads to an increase in the specific interfacial area at a given saturation. The increase in the specific interfacial area becomes more prominent around the maxima where the difference between the interfacial area for reconstructions with porosity of 0.6 is nearly 30% larger than the interfacial area for the reconstructions with porosity of 0.3.
It must be noted that the liquid-vapor interfacial area would depend on the wettability of the pores and the pore sizes. Therefore, this property would strongly depend on the type of CL used and the data should be treated with care when trying to extrapolate to different CLs.

Conclusions
Numerical simulations on CL stochastic reconstructions, generated using an overlapping sphere algorithm, were used to determine a general equation to estimate dry and wet effective diffusivity in CLs. First, a representative CL microstructure was obtained by comparing statistical correlation functions for the stochastic reconstructions to those obtained from the FIBSEM reconstruction of a thin, low loading CL. It was found that stochastic reconstructions with r d = 40 nm had two-point correlation function, chord length function, pore size distribution, specific interfacial area and mean chord length similar to those from the FIBSEM reconstruction thereby, warranting the use of these reconstructions to study the gas diffusivity in the CL. The detailed comparison of the statistical functions also provides a method to determine the representativeness of the stochastic reconstructions for fuel cell porous media, which has largely been ignored in prior numerical studies.
Then, gas diffusion was studied in the stochastic CL reconstructions with varying porosities under dry conditions. Analysis of the percolating pore network at different porosities showed that for porosities below 6% there was a loss of percolating gas network in the CL reconstructions. Effective diffusivity was isotropic for the stochastic reconstructions and increased with an increase in the porosity. A correlation for the effective diffusivity of the CL and its porosity was developed based on percolation theory. An empirical correlation was developed for the Knudsen radius based on the analysis of the pore size distributions for the stochastic reconstructions and used to compute the bulk diffusion coefficient in the pores which was used in the correlation. The percolation threshold was computed from the analysis of the percolating pore network and found to be 0.05 while the exponent for the correlation was computed to be 1.90. The effective diffusivity predictions for the stochastic reconstructions under dry conditions were also compared to experimental and numerical values from literature and found to be in good agreement.
To study the impact of local saturation on the gas transport in the CL, a novel nucleation based water intrusion algorithm is presented. A cluster based full morphology algorithm is used to track the liquid water interface. The liquid-vapor interfacial area which is critical to estimate the evaporation rates in macro-homogeneous models was calculated as a function of saturation for the different CL reconstructions. The liquid-vapor interfacial area increased to a maximum up to a saturation of nearly 75% and then decreased to zero at 100% saturation. Gas diffusion simulations were carried out on the partially saturated CL reconstructions with dry porosities in the range of 0.3-0.6. To account for the change in the pore morphology and Knudsen effects, pore sizes were recalculated assuming liquid water to be impervious to gas. The predicted effective diffusivity decreases with an increase in saturation for all the reconstructions. At a given saturation, the ratio of wet effective diffusivity to dry effective diffusivity were similar for all reconstructions irrespective of their porosity. The predicted wet effective diffusivities were used to develop a correlation valid for both dry and wet effective diffusivity in the CLs.